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I. INTRODUCTION AND SUMMARY 


The following summarize the research projects undertaken on NASA 
Grant NGR 05-017-031. Detailed research results are presented in Sec- 
tions II through V. 

PROJECT': Application of Vector Magnetometer to Aircraft Control 

and Navigation . " ~ 

INVESTIGATORS: W. R. Dunn and R. Pietila 

Recent technological improvements and corresponding cost reductions 
in 3-axis (i.e., vector) magnetometry led to investigation by Ames and 
University of Santa Clara researchers into use of this device in avionics 
flight control and navigation applications. Early analytical work showed 
that vector magnetometer outputs could be used to determine aircraft 
change provided that one attitude angle was known. These results were 
confirmed by successful correlation of attitude computed from magnetometer 
data taken during a NASA Convair 880 test with actual aircraft attitude. 

To demonstrate the viability and uncover the limitations of the theory, 
development of a solid state remote magnetic indicator was undertaken. 

The effort involved development of a microprocessor based system employing 
a 3-axis magnetometer and evaluation of the system using rest apparatus 
constructed ar the University. The background and results of this work 
are described in detail in Section II of this report. Results of the re- 
search were also published in [1] and [2] . 

PROJECT: Application of Earth Electric Field Phenomena to Aircraft 

Control and Navigation. 

INVESTIGATORS : W. R. Dunn, C. Keller, W. Keller, R. Wilson, R. Yarbrough 

Hill's [3] late 1972 paper on use of the earth's electric field in 
flight control prompted NASA and University researchers to undertake in- 
vestigation of earth field applications to large aircraft control and navi- 
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gation. 

The principle underlying Hill ' s work is the fact that in clear weather, 
the earth's surface terminates a large (approximately 300v/m at the surface) 
electric field which is close to the local vertical. Hill then successfully 
demonstrated that airborne (RPV) measurements of this field could be used to 
provide pitch/roll control. 

University investigations were undertaken to -evaluate the use of this 
phenomena 'in bad (e.g. IFR) weather conditions. The principal finding of 
this investigation was that this field (direction) can be significantly 
altered in and near cloud structures particularly the cumulonimbus cloud 
structure. Section III provides a review of the basic earth electric field 
phenomena and develops the basic field model for the cumulonimbus structure. 
The basic conclusion of this work is that the variability of the effects of 
the cumulonimbus on the clear weather electric field made navigational use 
of this field in all flight conditions a highly unreliable proposition. It 
is further proposed that electric field measurements might possibly better 
serve in a low cost system for identifying the direction and change intensity 
of cumulonimbus clouds during flight. 

PROJECT: Electronic Aircraft Cabin Noise Suppression 

INVESTIGATORS: W. R. Dunn and W. Keller 

Section IV describes an extension of Olson's [4] pioneering work m the 
use of active compensation as a means of reducing acoustic noise. This ex- 
tended work provides a theoretical basis on which wide band acoustic noise 
can be significantly reduced in aircraft cabins. 

PROJECT: SIRU (Strapdown Inertial Reference Unit) Technical Support . 

INVESTIGATORS: W. R. Dunn, W. Keller, C. Keller 


University- researchers worked with government and (government) contractor 
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personnel in the installation and checkout of the Carco motion simulator 
facility developed for use in the SIRJ program. This effort consisted of 
Carco table orthogonality measurements, adjustment and calibration of re- 
cording electronics, and specification of Carco Table/ Sigma 7 computer in- 
terface. 
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Aircraft Attitude Measuremen 
Using a Vector Magnetometer 



BLANK. NOT 


preceding page 

tabl¥^^pM« BLank Noi 


Acknowledgment iii 

Abstract ^ iv 

I A VECTOR AUTOPILOT SYSTEM 

1-1 Introduction 1 

1-2 Attitude Determination ... 2 

1-3 Attitude Determination Employing Magnetic 

Field 'Components 6 

1-4 A Possible System Configuration 8 

1-5 Other Considerations 15 

1- 6 Conclusions 19 

II AN ATTITUDE, INDEPENDENT REMOTE MAGNETIC INDICATOR 

2- 1 Introduction 2 0 

2-2 An Algorithm to Compute Aircraft Heading ... 21 

2-3 Mechanization of the Heading Algorithm .... 27 

2- 4 Conclusions 28 

I’ll DESIGN OF A MICROPROCESSOR BASED HEADING INSTRUMENT 

3- 1 Introduction 30 

'3-2 Hardware Design Considerations 31 

3-3 Software Design Considerations 37 

3-4 Design of Subroutines 39 

3- 5 Conclusions 50 

IV HEADING INSTRUMENT ERROR ANALYSIS 

4- 1 Introduction 54 

4-2 Sensor Errors 55 

4-3 Analog Subsystem Error Analysis 74 

4-4 Processing Errors 90 

4-5 Measurement Error •Summary 97 

4-6 Sample Error Analysis 99 

4-7 Conclusions j_02 


vi 



Table of .Contents (Continued) 

V LABORATORY EVALUATION OF THE ATTITUDE INDEPENDENT 
REMOTE MAGNETIC INDICATOR AND HEADING INSTRUMENT 

5-1 Introduction 104 

5-2 Test Apparatus 105 

5-3 Heading Measurements With No Offset Correction 109 

5-4 Heading Measurements to investigate Orthogon- 
ality Error 113 

5-5 Conclusions 114 

Appendix A: Instruction Set of the Signetics 2650 

microprocessor chip 127 

Appendix B: Assembly Language Program 130 

Appendix C: Table Generating Programs 156 


vii 



CHAPTER I 


A VECTOR AUTOPILOT SYSTEM 


1-1 INTRODUCTION 


An essential requirement of an aircraft attitude con- 
trol system is that deviation of the body axes relative to 
a reference axes frame must be sensed. In addition, to 
overcome the ever-present possibility of errors or failure 
of the sensors, various configurations of redundant sensors 
are usually employed to assist in detection and correction 
of errors. To this end, there has been a continuing effort 
to improve existing sensors, to develop ’new sensor config- 
urations , and to develop new sensor devices . 

This chapter discusses the role of a vector magneto- 
meter 1 as a new instrument for, aircraft attitude determin- 
ation. Although magnetometers have played a role m the 
attitude measurement of missiles and satellites [Ref. 1-1] , 
there is an apparent lack of application in aircraft systems 
By providing independent measures of attitude, the solid 
state vector magetometer sensor system can not only assist 
m improving accuracy and reliability of existing systems 
but can also reduce component count with obvious benefits 
m weight and cost. Additionally, since a large number of 
aircraft heading reference systems depend on measurement 
of the Earth’s magnetic field, it can be shown that by sub- 
stituting a three-axis magnetometer for the remote sensing 
unit; both heading and attitude measurement functions can 
be derived using common elements, thereby further reducing 
the component count. 


Aviation use to date has been essentially scalar magneto- 
metry . 


1 
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To investigate the feasibility of the above system, this 
chapter will proceed by developing a technique to determine 
attitude given magnetic field components. Sample calcula- 
tions are then made using the Earth’s magnetic field data 
acquired during actual flight conditions. Results of these 
calculations are compared graphically with measured attitude 
data acquired simultaneously with the magnetic data. The 
role and possible implementation of various reference angles 
are discussed along with other pertinent considerations . 
Finally, it is concluded that the Earth's magnetic field as 
measured by modern vector magnetometers can play a signifi- 
cant role in attitude control systems. 

1-2 ATTITUDE DETERMINATION 


Coordinate systems are usually defined by orthogonal 
right-handed sets of three unit vectors. An example of such 
a set is illustrated in Fig. 1-1 where the orientation of 
the body fixed frame used in this paper is delineated. An- 
gular- rotations are conventionally defined as rotations in 
the plane normal to a unit vector with the positive sense of 
rotation defined by the right-hand rule [Ref. 1-2] . 

To derive relationships of attitude variations as a 
function of magnetic vector component variation , we can pro- 
ceed by considering matrix representations of an orthogonal 
transformation. If Hx, Hy, and Hz are the magnetic compon- 
ents measured at a desired airframe attitude and Hx’, Hy', 

and Hz' are the components measured after any rotation of 

T 

the body, vector H' = [Hx 1 Hy' Hz'] can be related to 

T 

vector H = [Hx Hy Hz] by an orthogonal linear transform- 
ation H' = AH. Here A must satisfy the orthogonality con- 
. . T T 

dition AA = I, where A is the transpose of A; addition- 
ally, the determinant of A must be unity [Ref. 1-3, 1-4]. 



Fig. 1-1 AXIS ORIENTATION 


Rotations about the z axis in Fig. 1-1 result m yaw 
deviations (iji) and in new components (H’) , as shown by 


Hx’ 


cos f sm i|j 0 


HX 

Hy' 

Hz' 


-sin cos i|> 0 

0 0 1 

! 

Hy 

Hz 


Similarly, independent rotations about the y axis and 
the x axis result in pitch (e) and roll ($) dependent varia 
tions in the measured H components, as shown by 
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The effect of a combined rotation can be expressed by 
using the product of the transformation matrices. In add- 
ition, if the rotations are small, the total rotation exper- 
ienced by applying sequential rotations is independent of 
the order in which the rotations are performed (Ref. l-3,l-4]. 
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Assume that the angular variations 0, if», and <j> are 
small enough so that the small angle approximations 

sin 8 = 0, sm ip - ip , sm — <f> , 

cos 0 = cos Ip = cos <{> = 1 

can be made. Then, if the .products of small angles (in rad- 
ians) can be assumed to be much smaller than the angles 
alone, the expression reduces to 

"Hx'l 1 ip — 0 ~1 THx" 

Hy' = -<p 1 <p Hy 

_ Hz 1 J 9 -<f. 1J [hzJ (1-5) 

Further modifications m the form of the matrices result m 

Hx* [~-Hz Hy 0 "1 0 1 [Hx~ 

Hy 1 = 0 -Hx Hz ip t Hy 

_ Hz 1 J Hx 0 -HyJ L <f> J LHz J (1-6) 

By subtracting, we arrive at an expression for the dif- 
ference in H components as functions of angular deviation. 

'Hx'l f Hz 1 AHx "I f-Hz Hy 0] re“ 

Hy' - Hy = AHy - 0 -Hx Hz 

.Hz 1 j L Hz J |_ AHz J Hx 0 -HyJ L ^ J (1-7) 

It is significant to note at this point that the trans- 
formation matrix is singular implying than solutions for 
0, i|», and'<j> are not independently available. 
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1-3 ATTITUDE DETERMINATION EMPLOYING MAGNETIC 

FIELD COMPONENTS 

A given . orthogonal set of three unit vectors can be dis- 
placed m Euclidean space by rotating the system through any 
angle 6 about a directed rotation axis. It is also customary 
to represent this rotation vectorially as a directed line 
segment whose length is proportional to the rotation angle. 
This rotation is analogous to the rotation experienced by the 
body fixed frame of Fig. 1-1 as the aircraft experiences 
combined pitch, yaw, and roll variation. During flight the 
body fixed set rotates about this rotation axis assuming new 
(possibly erroneous) attitudes in space. The task of the 
attitude sensing system is to provide measures of compounded 
pitch, yaw, and roll that would result m the same attitude 
assuming that the rotations occurred sequentially about the 
x, y and z axes rather than the actual rotation axis . 

It was shown in the previous section that a compounded 
rotanion of an orthogonal set can be described by a product 
of respective . transformation matrices. Additionally ir was 
noted that for small angular rotations the order of multi- 
plication is unimportant. Using the relationships of (1-7) , 
expressions for the angular deviations m terms of measured 
magnetic vector components can be derived. 

AHx = -Hz9 + Hy<|> (l-8a) 


yields 

9 = (Hyi[> - AHx)/Hz (l-8b) 

* = 


( AHx + Hz0)/Hy 


(l-8c) 
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Similarly, 


AHy = -Hxi{i + Hz<j) 


(l-9a) 


yields 


* = - AHy) /Hx 

9 = (AHy + Hx^)/Hz 


(l-9b) 

(l-9c) 


and 


AHz = Hx0 - Hy<f> 


(l-10a) 


yields 


9 = (AHz + Hy£)/Hx 
<i> = (Hx0 - AHz) /Hy 


(l-10b) 

(l-10c) 


Assuming that Hx, Hy and Hz are nominal vector compon- 
ents as measured in a reference attitude and that Hx' , Hy 1 
and Hz' are new field components at the new attitude, then 
AHx - Hx' - Hx, AHy = Hy' - Hy, AHz = Hz ' - Hz are expres- 
sions of the incremental changes m field components. Add- 
itionally, before using (1-8) , (1-9) or (1-10) to solve for 

attitude variations (pitch, yaw, or roll) , one additional 
angle from an auxiliary sensor 2 must be supplied. Using one 
additional angle of rotation (about any one axis) the remain- 
ing two rotations can then be- calculated. 


To illustrate this point, flight data acquired during 
the flight of a NASA flown Convair 900 instrumented with a 
three-axis magnetometer and a Litton inertial navigation 
system were used to calculate roll, pitch, and yaw. 


2 It was noted following (1-7) that a unique solution for 
attitude variation is not possible using magnetic field 
data alone. * 
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Attitude variation about each of the three axes was calcul- 
ated using measured magnetic field components supported by 
one angle from the inertial system. The results of these 
calculations are plotted in Figs. 1-2 through 1-4. 

It is significant to note that the rotations shown 
occurred simultaneously (i.e., time base is the same for all 
three figures) . The flight was at an altitude of approxi- 
mately 5000 ft at an airspeed of approximately 250 nmi/h. 

Although the data used to plot the attitudes shown in 
Figs. 1-2 through 1-4 were not acquired specifically for this 
purpose, the correlations in measured and calculated atti- 
tude clearly show that, within the limits of instrument: ac- 
curacy, signals proportional to attitude variation can be 
derived using flight data. 

1-4 A POSSIBLE SYSTEM CONFIGURATION 


Since the intent of this chapter is to introduce the 
notion that magnetometer technology has advanced to the point 
where three- axis magnetometers can be incorporated in air- 
craft attitude sensing systems on a cost effective basis, the 
system discussion will be limited in scope no describing a 
possible combined heading and attitude measurement method. 

Heading references fall into three classes; 1) those 
that depend on rhe Earth's magnetic field, 2) those that de- 
pend on the use of low-drift gyroscope to retain a preset 
azimuth, and 3) those (gyrocompasses) that depend on sensing 

t 

the Earth's rotation [Ref. 1-5). By far the greatest number of 
aircraft heading systems depend on the Earth's magnetic 
field, although many of these include gyroscopes to' improve 
the performance characteristics. 



ROLL { DSGRE 





PITCH CALCULATED USING MAGNETIC FIELD DATA 


PITCH MEASURED 








0 24 48 72 96 120 144 168 192 216 240 

TIME (SEC) 


EICURE 1-4 YAW AXIS 




12 


A popular system combination (with no gyro) is to 
combine a pendulous remote magnetic sensor and a synchro 
receiver m a null seeking circuit. The philosophy being to 
attempt to measure only the horizontal component of the 
Earth's magnetic field and to swing the receiver into align- 
ment with it. Under acceleration, departures of the sensor 
unit from the horizontal result in angular heading errors e 
[Ref . 1-5] . 

e = (aH/g) tany sine 

where aH is the horizontal acceleration, g is the accelera- 
tion due to gravity, 0 is the angle between the acceleration 
vector and magnetic north, and y is the magnetic field dip 
angle; arctan (vertical f leld/horizontal field) . 

Accuracy of this system can be improved by incorporat- 
ing a strapped-down solid state magnetic sensing unit (free 
of acceleration errors) that measures and displays the angle 
of the Earth's horizontal magnetic component relative to the 
aircraft. This system can be implemented as follows: 

1) Determine the direction of the magnetic vector F rela- 
tive to the sensors (and the airframe) , by measuring the x, 
y and z components (Figs. 1-1 and 1-5) . The direction co- 
sines cosa, cos 0 , cosy are the cosines of the angles a, 0, y 
between the magnetic vector and the positive x, y and z 
axes. Additionally, 

cosa = x/(x 1 2 + y 2 + z 2 ) 2 

cosg = y/(x 2 + y 2 + z 2 ) 2 

cosy = z/ (x 2 + y 2 + z 2 )* 
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2) Using either a vertical reference 3 or knowledge of 
aircraft attitude, we -can effectively rotate the body axes 
such that the x-y plane is horizontal (see Chapter II) . 

3) Simple application of direction cosines will yield the 
direction of magnetic north in the aircraft's x-y plane. 


X 


z 


Fig. 1-5 FIELD VECTORS AND DIRECTION COSINES 

Although the preceding discussion implies that heading 
can be determined by using a strapped-down magnetometer, 
there remains the problem of attitude determination. An- 
other widely used system for obtaining a heading reference 


3 Not necessarily derived inertially (Ref. 1-11]. 
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Is to combine the relatively excellent short term stability 
of a directional gyroscope with the long term stability of 
magnetic field measurements. By slaving the directional gyro- 
scope to the magnetic field [Ref. 1-5, sec. 10.4.7] , gyroscopes 
with relatively large free drift error can be used s to pro- 
vide an excellent heading reference. 

Replacement of the pendulous remote sensing unit of this 
type of system with a strapped-down vector magnetometer would 
result in both heading and attitude information on a contin- 
uous basis. This combination would operate as follows: 

1) The system is initialized by determining a reference 
attitude (perhaps by using a primary inertial attitude sys- 
tem) . 

2) The angular position of the horizontal magnetic field 
component is computed as above and used to slave the direct- 
ional gyroscope. 

3) The directional gyroscope, with relatively good short 
term stability (devices with free drift of less than 0.5 
deg/h have been designed) , is used to determine yaw (ij>) 
errors . 

4) For small angle deviations, (1-8) , and (1-9) , and 
(1-10) can be employed to recalculate aircraft attitude. 

The process loops back to step 2) closing the loop on a 
combined attitude and heading reference system. 

The sampling frequency required to maintain an accept- 
able level of error is of course determined by the aircraft 
performance expected (angular rates) and by the gyro error 
(drift rate plus errors due to additional sources such as 
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gyroscope tilt from vertical) . The overall system is such 
that heading can be determined as before with errors due to 
sensor departures from horizontal substituted for long term 
accumulation of attitude uncertainty {this can be corrected 
by looping to step 1) at a frequency dependent on error 
rates) . Additionally one gains measurements of attitude 
with minimal computation and replacement of a mechanical re- 
mote sensing unit with a solid state strapped-down magneto- 
meter sensor. 

1-5 OTHER CONSIDERATIONS 

The characteristics of the Earth's magnetic field and its 
variations have long been established [Ref. 1-6-1-10] . Since 
the field is to be used as a reference in the attitude 
measurement scheme, there is a need here to discuss its ad- 
verse characteristics. Although the field does experience 
variation, most of the variation is either in amplitude 
(ionospheric contributions) or has time constants that make 
the variation negligible (secular variation) . 

In traversing local anomalies, there will, however, be 
deflections m the ambient field due to the additive effect 
of local dipoles or monopoles. The effect of local terrain 
caused anomalies can be visualized by picturing the main 
field vector oriented m space with a second modulating 
vector rotating at its tip. Maximum angular error would 
occur when this modulating vector has maximum magnitude and 
is positioned at right angles to the mam vector. 

To illustrate the effect of local anomalies one can cal- 
culate the level of anomaly required to cause an error. 

Since the Earth's mam field is typically in the order of 
0.50 G it is readily apparent that a local anomaly of ap- 
proximately 0.01 G at right angles to the local field is 
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required to cause an error of 1 deg. Furthermore, the 
local anomaly would have to be aligned with one of the air- 
craft body axes to result in one degree of attitude error in 
any one axis. Fortunately, anomalies with components of 
this magnitude positioned at right angles to the main field 
are extremely rare. In addition, the anomalies are local- 
ized over ore bodies or other geophysical irregularities, 
have magnitudes that diminish as the cube of altitude , and 
tend to average to zero over relatively short distances. In 
summary, the probability of encountering an anomaly that 
would cause as much as a 1 degree error is relatively small. 
The error, if introduced, will be short lived and, unlike 
drift error, will average to zero. 

Fundamental to a magnetic field referenced system is 
the ability to measure orthogonal components of the field 
vector. Precision and accuracy of measurement of the com- 
ponents is of course specified by the desired control spec- 
ifications . 

Since the Earth's magnetic field varies m magnitude on 
a global basis between 0.3 G and 0.6 G (30,000 gamma to 
60,000 gamma), it is apparent that full scale measurements 
of 0.6 G can be expected. Sensors mounted at right angles 
to the field will monitor no measureable field and thus de- 
fine the lower limit of measurement to be zero. For the con- 
tinental United States the declination varies between 60 and 
80 deg, resulting in a range in horizontal component of 0.15 
to 0.25 G with vertical component in the range of 0.4 to 0.55 
G. Heading variations (yaw) result m changes of the hori- 
zontally sensed field components and would specify the maxi- 
mum precision required. In addition, flight at 45 deg ± (n 
x 90 deg) (where n is any whole number) with respect to mag- 
netic north results in minimum sensitivity of the x and y 
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axes measurements. In this case sensor inputs would range 
between 0.106 and 0.177 G with minimum field at the north. 
Assuming the preceding ambient measurements , variations in 
component magnitude of approximately 0.0180 to 0.0305 G/deg 
for small angle variations can be expected. 

A brief survey of commercial magnetometer manufacturers 
reveals that triaxial magnetometers that measure from zero 
to 0.6 G with linearities of 0.5 percent, noise less than ±1 
mG and sensitivies of at least 2.5 V per 600 mG are current- 
ly available. In addition, these devices have a bandwidth 
of direct current to at least 500 Hz and are rated to have 
less than 1 deg error in orthogonality. 

From a precision standpoint, it is apparent that 'vari- 
ations in yaw for this worst case situation can be sensed ro 
better than 0 . 1 deg with currently available magnetometer 
technology. The sensor technology required to implement an 
attitude sensing system of reasonable specifications is avail- 
able ‘(more detailed analysis is presented m Chapter III). 

Although the preceding calculations indicate that for 
small angular variations attitude can be calculated using 
measured magnetic data, there is a need to consider the 
effects of larger finite rotations. In this case the small 
angle assumptions would not be valid and an Euler transform- 
ation would have to be made. Measurement of three axes of 
field components could be used to develop the direction co- 
sines required to determine the orientation of the axis of 
rotation, the angular rotation about it, and the three angu- 
lar rotations of pitch, roll, and yaw. 

For the special case where the axis of rotation aligns 
with the magnetic vector, there would of course be no 
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measured component changes . ** By measuring the attitude of a 
second vector (not in alignment with the magnetic vector) , 
we could resolve the ambiguous situation cited above and pro- 
vide additional redundancy. 

The optimum auxiliary vector would be one that could be 
sensed without using inertial devices. The Earth's electric 
field can be considered. The mam reason for considering 
this field as a means of providing an auxiliary angular ref- 
erence is that the resultant system has the potential of be- 
ing completely solid state. The electric field vector can 
be used to determine attitude variation m a manner analogous 
to the magnetic vector system. Inherent limitations of each 
single vector system can be obviated if the vectors are not 
coincident. 

Although Hill [Ref. 1-11] reported success in controlling 
pitch and roll using the electrostatic field alone, comments 
by MarJcson [Ref. 1-12] indicate that the electrostatic field 
is not always a reliable vertical reference. Employment of 
the electrostatic field for this attitude measurement system 
is limited to augmenting the magnetic field measurements by 
eliminating ambiguity of motion around the magnetic vector . 

The requirement of vertical electrostatic field is thus remov- 
ed and replaced by a requirement that the field direction is 
relatively stable. 

By using two independently derived vectors we have suf- 
ficient data to obviate the ambiguity just cited and we have 
the potential of providing redundancy as well. 


**An example of this would be yaw rotation while flying 
straight and level over the magnetic poles or roll rot- 
ation while flying towards a pole at the magnetic equator. 
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1-6 CONCLUSION 

This chapter has identified a novel method of measuring 
aircraft attitude using relatively inexpensive, well devel- 
oped instrumentation. It has recognized that magnetic field 
sensing systems have been used to some extent in attitude 
sensing and control of space vehicles; it has also suggested, 
however, that with appropriate support, magnetometers can 
find increased application m aircraft attitude measurement 
systems . 

This claim is corroborated by actual flight test data. 
Magnetometers have evolved to a point where three axis mea- 
surements of the Earth's magnetic field can be made with 
sufficient precision and accuracy to enable measurement of 
small angle attitude variations . 

This chapter has also discussed a possible system con- 
figuration combining heading determination and attitude 
measurement functions. By replacing the conventional re- 
mote sensing unit with a three- axis magnetometer, it has 
been suggested that both functions can be obtained with the 
hardware required previously for heading measurement alone. 

As with any system, there are limitations imposed. The 
main limitation for a vector magnetometer system seems to be 
the inability to sense rotations around the magnetic vector 
itself. This problem is not unlike the ambiguity experienc- 
ed by magnetic heading systems at high latitudes. By judic- 
iously incorporating auxiliary instruments, not only can 
the ambiguities be removed but a degree of redundancy can be 
added while still maintaining a cost and weight advantage 
over comparable systems. 



CHAPTER II 


AN ATTITUDE INDEPENDENT REMOTE MAGNETIC INDICATOR 
2-1 INTRODUCTION 

Preliminary investigation [Ref. 2-1] revealed that air- 
craft attitude' can be calculated using measurements of earth's 
magnetic field vector and a single auxiliary rotation angle. 

An algorithm to compute the two remaining aircraft rotational 
angles was developed. Using flight data, it was demonstrated 
that an excellent correlation m computed versus actual air- 
craft attitude could be achieved. In addition to providing 
measurements of the magnetic field for redundant attitude* 
computations (to improve accuracy and reliability of existing 
autopilot systems) , it was noted that the vector magnetometer 
could substitute for the remote magnetic sensing unit. In. 
this manner both heading and attitude measurements could be 
derived using common elements with obvious benefits in weight 
and cost. 

This chapter discusses the mechanization of a micropro- 
cessor based computer system that uses a three axis magneto- 
meter plus gyro data to compute heading. The magnetometer is 
a three axis solid state device that can be mounted in a 
strapped down configuration resulting in an attitude independ- 
ent remote magnetic indicator. Gyro measurements of pitch and 
roll angle plus three axis magnetic measurements are used by 
the algorithm to compute aircraft heading. The system can 
function independently to compute heading or by simply in- 
creasing the stored program' could implement the attitude com- 
puting algorithm of [Ref. 2-1] as well. 

The chapter proceeds by developing an algorithm to com- 
pute aircraft heading using the strapped down magnetometer 
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and two gyro measured angles. Practical aspects of designing 
the system including, both hardware and software are then pre- 
sented. In addition, the limitations m instrument accuracy 
and operation as determined by sensor errors, signal process- 
ing errors, arithmetic precision and computation speed are 
discussed. Considerable computational capability inherent in 
the system enables minimization of systematic errors. It is 
demonstrated that inexpensive sensors can be employed with 
offset and orthogonality errors compensated by microprocessor 
programming. Finally, it is concluded that a microprocessor 
based computer with a solid state magnetometer can play a sig- 
nificant role m aircraft instrumentation. 

2-2 AN ALGORITHM TO COMPUTE AIRCRAFT HEADING 

Coordinate frames are usually defined by orthogonal 
right-hand sets of three unit vectors. An example of such a 
set is illustrated in Fig. 1-1 where the orientation of the 
body fixed frame used in this chapter is delineated. The ref- 
erence coordinate frame referred to in this chaprer is orient- 
ed with axes x and y in the horizontal plane and axis z verti- 
cal (z down is positive) . Pitch attitude angle (0) of an air- 
craft is defined [Ref. 2-2] as the angle between some prefer- 
red longitudinal axis and the horizontal reference. In this 
chapter, pitch angle is the angle between the x axis of the 
aircraft and the x-y plane of the reference axis set. Since 
angular rotations are conventionally defined as rotations m 
the plane normal to a unit vector with the positive sense of 
rotation defined by the right-hand rule [Ref. 1-2], we will 
define positive pitch angle (0) as the "nose up” or positive 
rotation about the y axis when the v axis is horizontal. The 
roll and yaw angles (4> and ¥) will then simply be rotations 
about the x and y axes respectively. 
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By aligning the three magnetometer axes with the respec- 
tive x, y and z axes of the aircraft, -we can measure magnetic 
field components of the aircraft at any attitude. For the 
trivial case where pitch (0) and roll (*) are both zero de- 
grees, Hx and Hy are the horizontal field components and we 
can compute yaw from the horizontal vectors as follows : 

H'i = cos 1 (Hx/(Hx 2 + Hv 2 ) 5 (2 -la) 


or 


= sin 1 (Hy/(Hx 2 + Hy 2 )^ 


(2-lb) 


We select either (2-la) or (2-lb) based on the relative mag- 
nitudes of Hx and Hy. By minimizing the numerator of the 
argument we guarantee that the inverse trigonometric operation 
results m an angle between zero and forty-five degrees with 
maximum sensitivity ensured. Heading is then computed using 
the signs of Hx and Hy to select the appropriate equation from 
Table 2-1. 



Table 2-1. Formulae to Compute Heading 

For most cases, the pitch and roll angles are not zero 
and inverse rotations are required to determine the actual 
horizontal field components Hx and Hy. Since any aircraft 
attitude can be represented as a sequence of rotations about 
each axis beginning at some reference attitude, we can 
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determine the reference Hx and Hy field components by perform- 
ing an inverse roll followed by an inverse pitch computation 1 . 


The inverse roll computation can be developed by consid- 
ering vector components of an arbitrary vector H in Fig. 2-1. 
The first set (x 2 , y 2 , z 2 ) represents the vector components 
measured in a reference orientation. The second set has com- 
mon origin and aligns with common x axis component. It is 
rotated (rolled) about the x axis resulting in new y and z 
values. We can describe vector H m both coordinate frames as 


H = x 2 • i 2 

+ y 2 

* 32 + 

z 2 ♦ k 2 

(2-2) 

A 

H = x 3 » i 3 

+ Y 3 

♦ 3 3 + 

Z 3 « k 3 

(2-3) 


Since the vector H is unique, we note that equations (2-2) and 
(2-3) are equal. Furthermore if we form dot products we solve 
for the horizontal components x 2 , y 2 , and z 2 in terms of 
the rotated values and the roll angle ($) . 

, From (2-2) we obtain 


H • 1 2 = x 2 (l 2 • 3-2) + y 2 (32 * i-2) + 2 2( k 2 


H • 


3-2 


= X 2 


(2-4a) 

(2-4b) 


and from (2-3) we obtain 


<N A A A A 

H - i 2 = X 3 U 3 * i 2 ) + Y3 ( j 3 • i 2 ) + Z'3( k 3 


H * i 2 = x 3 


(2-5a) 

(2-5b) 


1 Since pitch is defined as the angle between the x axis and 
the horizontal plane we can assume that at any heading, air- 
craft attitude results due to a pitch followed by a roll. 
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then 

Similarly, 


x 2 = x 3 


( 2 - 6 ) 


H * j 2 = Y 2 -X 3 (i 3 * j 2 ) + 73^3 * 3 2 ^ + z 3^ 3 


y 2 — y 3 cos $ - 2 3 sin <i> 


and 

/\ a ^ ^ ^ * 

H • k 2 - Z g — X 3 ( 1 3 * k 2 ) + Y3 (3 3 * k 2 ) + 23(^3 
z 2 = Y3 sm 4 + z 3 cos <p 
These expressions can be summarized as 


/S 

32 ) 

( 2 - 7 a) 
( 2 - 7 b) 

k 2 ) 

(2-8a) 

(2-8b) 


*2' 


"1 

,0 

0 



~x 3 -' 

y 2 

= 

0 

cos <i> 

-sin 

$ 

• 

y 3 

z 2 _ 


_0 

sm § 

cos 

4 > _ 


_ Z '3 . 



Fig- 2-1 AXES ROTATED IN ROLL 
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Similarly, considering an axis set rotated in pitch as shown 
in Fig. 2-2, we can express the reference set x^ , yj , z x m 
terms of the rotated set X 2 / z 2 as follows 


■*r 


cos 0 0 sin © 1 


~x 2 - 

y i 


O 
i — I 

o 

* 

y 2 

.21 . 


-sin 0 0- cos 0 J 


_ z 2 _ 



Fig. 2-2 AXES ROTATED IN PITCH 

Finally, if we assume that the axis set subscripted with 3 
represents components of Earth's magnetic vector measured at 
an arbitrary aircraft attitude, we can derive the magnetic 
components (Hxh, Hyh, Hzh) in the horizontal plane for a 
given heading 
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” Hxh' 


cos 0 0 sin 0 ' 


H 1 
O 

O 


“Hx 3 " 

Hyh 

=* 

0 10 


■ 0 cos $ -sin. $ 

♦ 

Hy 3 

_Hzh _ 


_-sin 0 0 cos 0 _ 


_0 sin © cos 


_Hz 3 . 


(2-lla) 


'Hxh' 


cos 0 (sin 0 sin ©) sin 0 cos 


'Hx 3 ' 

Hyh 

= 

0 cos © --sin -$ 

• 

Hy 3 

_Hzh_ 


_-sin 0 (cos G sin <f>) cos * cos Q_ 


_Hz 3 _ 


( 2-llb) 


The algorithm to be implemented with the microprocessor 
would therefore require operations as 'Outlined in Fig. 2-3. 
Details of programming method, modifications to the above 
equations to facilitate programming and computation speed 
versus accuracy tradeoffs are discussed in following sections. 


DO FOREVER 


Measure, Digitize and Store 
Hx , Hy , Hz , 0 and © 


Correct Sensor Errors 


Compute the Horizontal Field 
Components Using Equation 2-11 


Compute Heading Using Equation 2-1 


Display Heading 


Fig. 2-3 LOGICAL OPERATIONS REQUIRED TO COMPUTE HEADING 
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2-3 MECHANIZATION OF THE HEADING ALGORITHM 

A. General Considerations 

To evaluate the performance of an integrated system ex- 
perimentally, an instrument was designed to implement the al- 
gorithm developed above. Several approaches were considered 
to implement the heading instrument for experimentation: 

1) A minicomputer implementation incorporating an 
HP-2100 minicomputer supported by peripheral interface and an- 
alog circuitry. Programming of the HP-2100 would have enabled 
the computer to control multiplexing and processing of sensor 
data as suggested by Parish and Lee [Ref. 2-3] . 

2) A hybrid system composed of a remote data acquisition 
system to collect data from sensors for subsequent processing 
by a computer (possibly an HP-2100) . 

3) A digital/analog electronic implementation incorpor- 
ating the design of a special purpose computer to perform the 
required functions of a heading instrument. 

The first two approaches were abandoned since it was de- 
sirable to perform the experiments at various locations remote 
from a computer facility and to have data available immediately 
without having to rely on off-line computations at a later date. 
The design task then evolved to the design of a special pur- 
pose computer system to implement the algorithm, provide a 
means for evaluating the performance of the proposed algorithm 
and to allow modifications to the system if required. 

B. Design Criteria 


Having decided on the general approach to implementing 
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the algorithm it became -necessary to consider the performance 
criteria desired of the instrument. 

1) Accuracy 

As a design goal, an absolute accuracy of +1.0° in 
heading uncertainty was selected for the laboratory 
implementation. This accuracy is compatible with 
commercially available heading systems. 

2) Computation Speed 

The bandwidth of the system is determined mainly by 
the computation speed of the computer 2 . As a design 
goal, complete heading updates, once per second was 
established. 

3) Flexibility 

l 

A desirable feature of the laboratory evaluation 
instrument was considered to be flexibility. Re- 
visions or additions to the algorithm as predicted 
■ by experimental data should be incorporated with 
minimal redesign of the instrument. 

2-4 CONCLUSIONS 


An instrument designed to implement the heading algorithm 
developed above uses a three axis magnetometer to measure mag- 
netic field data m the vicinity of an aircraft. Since the 
magnetometer proposed is a solid state three axis fluxgate de- 
vice and is permanently mounted m a strapped down configura- 
tion, the implementation results m an attitude independent 


2 The response times of the various sensors and analog circuitry 
are orders of magnitude greater than the desired one second 
sample interval. 
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remote magnetic indicator 3 . 

Several factors will contribute to system inaccuracy. Al- 
though the major error sources can be evaluated mathematically 
j[ Chapter IV) , there is a need to evaluate the implementation 
experimentally. Systematic errors that arise can be reduced 
by instrument computation. This capability (inherent with a 
computer based system) enables incorporation of less expensive 
sensors in the heading instrument with less concern with fact- 
ors such as temperature regulation, sensor orthogonality and 
sensor offset 4 . 

Since the algorithm can be implemented using a micropro- 
cessor as the major computer element, the resulting instrument 
will have inherent computation capability, be small in size 
and consume relatively little power. These factors make the 
instrument an ideal device for aircraft application where the 
need for redundant distributed processing capability is in- 
valuable . 


3 Current remote magnetic indicators are pendulous and rely on 
gravity to enable measurements of the horizontal magnetic 
vector (not attitude independent) . 

4 Assuming that the sensors have repeatable or measurable ' char- 
acteristics, algorithms can be developed to correct previous- 
ly measured erroneous data. 
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DESIGN OF A MICROPROCESSOR BASED HEADING INSTRUMENT 
3-1 INTRODUCTION 

Progress in device and component technologies during the 
1970 ' s- has led to an assortment of sophisticated integrated 
circuits (IC) devices [Ref. 3-1] which enable the design of in- 
struments with a high degree of sophistication and accuracy. 

Of these devices, the microprocessor has to date been the most 
exploited component m industrial control and instrumentation 
applications [Ref. 3-2 through 3-7] . There have been many 
papers presented addressing the general application and feas- 
ibility of applying microcomputers to particular design tasks 
[Ref. 3-8 through 3-21] . 

Although much of the literature to date on microproces- 
sors has addressed the design of commercial products (usually 
rhe final result of a carefully orchestrated effort beginning 
with a market survey) , the design of a laboratory instrument 
for algorithm evaluation differs in design philosophy. In 
particular, the laboratory instrument is designed to evaluate 
a proposed algorithm under laboratory conditions. The tradi- 
tional benchmark evaluations and attempts to match the micro- 
processor to the application is not only difficult but unnec- 
essary. If the processor is much more powerful than necessary, 
the "overkill" is little noticed; but if an insufficiently en- 
dowed microprocessor is selected, the effects can be devastat- 
ing. Not only will the program be difficult to write and vor- 
acious of memory, it would be difficult to change to a more 
powerful microprocessor part way through the project. With 
these considerations in mind, a general purpose, flexible 
microprocessor with powerful architecture and instruction set 
the Signetics 2650 microprocessor [Ref. 3-22] was selected. 
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3-2 HARDWARE DESIGN CONSIDERATIONS 


‘The design of a microprocessor based system begins by 
considering the total system level block diagram to be imple- 
mented (Fig. 3-1). Inputs from five sensors including x, y 
and z axis magnetic data plus pitch and roll angles (Hx, Hy, 
Hz, 0 and *5) are to be multiplexed, sequentially sampled and 
converted to a digital representation prior to processing (ex- 
ecuting the algorithm developed above) . The main subsystem of 
Fig. 3-1, the central processing unit (CPU) , operates under 
control of instructions stored in the system memory and inter- 
faces with the input and output subsystems via data ports. 

At this early stage m the design, it is significant to 
note that the block diagram of Fig. 3-1 differs slightly from 
that of a classical discrete hardware solution. The input 
subsystem (composed of analog multiplexer, sample and hold, 
and analog to digital converter) differs from a conventional 
data acquisition in that it is devoid of a control section. 

The microprocessor will control the data acqui siti on sampling 
and conversion m addition to performing the arithmetic func- 
tion associated with the algorithm. 

Having established a tentative block diagram of the in- 
strument, the design continues by addressing relevant charac- 
teristics and limitations of each subsystem. These character- 
istics will then in term be considered m configuring the 
final system and program to be executed. 

1) The Analog Subsystem 

Composed of the analog multiplexer, sample/hold and an- 
alog to digital converter, the analog subsystem of Fig. 3-1 
affects both system accuracy and throughput rate. The 



THE ANALOG SUBSYSTEM 



Fig . 3-1. 


SYSTEM BLOCK DIAGRAM 


OJ 
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well-known Shannon theorem [Ref. 3-23, 3-24] on sampling 
theory defines one of the basic limits on throughput rate 
stating that the minimum frequency for sampling must be double 
the highest significant frequency of the signal, including the 
noise on the signal. This minimum frequency is necessary, the 
theorem states, if the sampled signal is to contain all of the 
information needed for undistorted reconstruction. At a lower 
sampling frequency ailiasxng can occur 1 . The minimum sampling 
rate for data to be used in this heading instrument (based on 
the design goal of Chapter II) then results in a system band- 
width of 30 hertz. The analog signals from each sensor are 
low pass filtered to reduce frequency content above 60 hertz. 

A survey of commercially available multiplexers, sample and 
hold modules and analog to digital convert modules (ADC) [Ref. 
3-25 to 3-28] reveals that subsystems with throughput -charac- 
teristics exceeding the requirements of a system sampled at 
one second intervals are readily available (pertinent specif- 
ications are discussed m more detail in Chapter IV) . The 
limiting parameter determining total system speed performance 
will .then be the execution time of the algorithm (a programm- 
mg consideration) . A further system consideration is the 
ability to adjust analog system offset and gain. These ad- 
justments are made using variable resistors (trim pots) con- 
nected to appropriate leads on the sample and hold and analog 
to digital converter modules . 

2) The Central Processing Unit (CPU) 

The central processing unit (Fig. 3-2) is composed of the 
microprocessor (Signetics 2650) supported by peripheral logic 
elements (Fig. 3-2) . Design of this subsystem involved medium 


^hac is, the sampled data derived from a sine wave of freq- 
uency f sampled at a frequency less than 2f can be fitted 
with sine waves of a frequency other than f. 



3 (8T‘J7 ) 2 ( 8T26) 
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and small scale integrated circuits using well-known [Ref. 

3-29 through 3-31] design techniques. To facilitate system 
development several features were included in the design of 
CPU subsystem (features that would not necessarily be re- 
quired in a production instrument). These include: 

a) System reset, single step and normal run mode oper- 
ation controlled by switches and logic elements . 

b) An RS-232 teletype interface is included to enable 
manual intervention and development capability 
during program development. The program was dev- 
eloped by loading and executing instructions into 
the random access memory (RAM) under control of the 
PIPBUG 2 program. 

3) The Memory Subsystem 

The memory subsystem (Fig. 3-3) was organized onto cards 
each with two thousand byte capability. In this manner system 
memory could easily be expanded (or reduced) in increments of 
2K bytes. The memory chips selected were organized as 256 
four bit words and feature pin for pin compatibility with 
commercially available random access (RAM) and programmable 
read only memory (PROM) chips. Program segments could then be 
developed in RAM and finally "burned" into PROM chips for a 
permanent, nonvolatile operation. In this manner the system 
development begins with IK bytes of memory devoted to the res- 
ident PIPBUG program (in ROM chips) with the remainder of mem- 
ory allocated as RAM for both program and scratch pad usage. 


2 S'ignetics tradename for the 2650 resident, loader and monitor 
program. 




Fig. 3-3. 
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As the program is developed, additional memory as added in 
increments of 2K bytes per card or 256 bytes on the card. 
Modifications to the program can be easily made using the 
PIPBUG program and teletype. 

4) The Output Subsystem 

For laboratory development the output subsystem of Fig. 

3-4 was designed to provide seven segment visual output of 
the aircraft heading with three significant digits displayed. 

To expediate the design cycle and to enhance system throughput 
rate, the outputs were designed as ports with latches and de- 
coder driver functions provided by hardware. In other applica- 
tions a hardware/software tradeoff could be made with the data 
decoding and driving implemented using table lookup and multi- 
plexing controlled by the CPU. 

3-3 SOFTWARE DESIGN CONSIDERATIONS 

■The general purpose processor selected to implement the 
CPU was designed to implement: programmed logic and to perform 
conventional computer operations. This heading instrument 
takes advantage of both areas. Since the instrument is act- 
ually a special purpose computer under control of a stored 
program, the functional specialization resides in the program 
rather than the hardware logic. Modifications can be made re- 
latively easily, satisfying the flexibility design goal of 
Chapter II. 

Having decided on the tentative hardware structure de- 
scribed in Section 3-2 above, the program development leading 
to the final listing in Appendix B proceeded as follows: 

1) Structured ' flow charts were developed depicting the total 

system operation as an ordered sequence of operations . Each 




Fig. 3-4 
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THE OUTPUT SUBSYSTEM 
















operation is identified as a separate subroutine which in 
turn can have "nested" subroutines of its own (Pig. 3-5) . 

2) System accuracy requirements were next investigated (dis- 
cussed in detail in Chapter IV) to ascertain the precision re- 
quirements 3 of the various subroutines. 

3) The respective subroutines outlined in 1) above were 
developed and implemented using a cross assembler program 
[Ref. 3-32] . Each subroutine was then loaded into the devel- 
opment hardware and "debugged" prior to total program inte- 
gration. The above program development depicts a top down 
strategy of program development [Ref, 3-33] and leads to an 
expedient system development with subroutines being individ- 
ually developed to yield a modular program construction. 

3-4 DESIGN OF SUBROUTINES 

The total program consists of an overall system program 
composed of nested subroutines. The discussion m this sec- 
tion is limited in scope to the design of the more-complex 
subroutines required to implement the solid state remote mag- 
netic heading algorithm. 

1) Subroutine "SAMP" (Fig. 3- 6a) 

The first portion of this subroutine is dedicated to the 

I 

control function of selecting an analog channel via the multi- 
plexer, sampling and holding the data, resetting and reading 
data from the analog to digital converter (ADC) . Prior to 
or during the programming of this section, data fields m 


3 This step is vital to determine whether the operations out- 
lined in 1) above are to be carried out in a single or multi- 
precision manner. 
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Power on reset of all registers 
and subsystems 


Compute the aircraft heading 


Fig. 3~5a. SYSTEM PROGRAM 




Sample all analog channels 

Compute the horizontal Hx field 

Compute the horizontal Hy field 

Compute the horizontal field 
vector 

Compute heading 
Output the data 


Fig. 3-5b. 


SUBROUTINE "MAIN 





















PRELOAD VARIABLES REQUIRED 
TO CORRECT FOR OFFSET ERROR 


SADD 


PERFORM A DOUBLE PRECISION ADD OR SUBTRACT 


RESTORE CORRECTED DATA TO TABLE 


RETURN 


Fig. 3-6b. SUBROUTINE "OFST" 




TEST 

FLAG 

ADD X. 

SUBTRACT 

DO DOUBLE PRECISION 


ADD ON SIGN 

DOUBLE PRECISION 

MAGNITUDE NUMBERS 

SUBTRACT 

RETURN 


Fig. 3-6c. 


SUBROUTINE "SADD 
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input ports 1 and 2 and output port 1 of Fig* 3-1 are allocat- 
ed. Control information ,is then passed to the peripheral mod- 
ule by writing control words to output port 1. Analog to dig- 
ital converter status and the 12 bit data field are sampled by 
reading input ports 1 and 2. 

Sensor outputs were biased at +2.5 Volts with transfer char- 
acteristics as depicted in Fig. 3-7a [Ref. 3-34]. The ADC 
selected for this laboratory instrument had a binary output 
data format related to analog input as shown m Fig. 3.7b 
[Ref. 3-35]. The second function of the sampling- subroutine 
"SAMP" was to convert data from a unipolar binary format to a 
sign magnitude format. Since the total transfer function from 
sensor input to ADC output (Fig. 3- 7a and b) indicates an off- 
set of 2.5 Volts or 1/2 the ADC output range, the sign magni- 
tude format can be generated as shown m Fig. 3-8. 



Fig. 3-7a. 


SENSOR TRANSFER CHARACTERISTIC 




Fig. 3-3 CONVERSION OF DATA 
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The third function of the "SAMP" subroutine was to reverse the 
sign of the Hx and Hy data (to correct a test fixture pro- 
blem) and to correct for sensor offsets. Although analog sub- 
system offsets are corrected by adjusting either the sample 
and hold module or the ADC, the independent sensors themselves 
have offsets 4 . Offset errors for the laboratory instrument 
were compensated by determining the offset correction term for 
each sensor (method described m detail in Chapter V) and then 
either adding or subtracting the term to the respective data 
during the sample subroutine. By characterizing the sensor 
errors 5 , actual datum could be improved further during this 
step . 


The final function of the "SAMP" subroutine was to cor- 
rect for sensor orthogonality error (subroutine "ORTH") . Al- 
though the sensors were physically aligned and specified to 
have orthogonality characteristic (Ref. 3-34] less than +1 
degree relative to the base coordinates , this nonorthogonality 
contributes appreciably to total system error (see error an- 
alysis m Chapter IV) . The physical misalignment of the sens- 
ors was determined experimentally (Chapter V) and determined 
to be mainly a misalignment of sensor x in the x-y plane as 
illustrated in Fig. 3-3. 

The actual data measured with the x axis sensor is then relat- 
ed to the true Hx and Hy values as 

Hx 1 = Hx Cos s - Hy Sin £ . 


4 With zero stimulus applied the sensors have a finite nonzero 
output. This error in the fluxgate magnetometer is a function 
of temperature, voltage and magnetic remanence in the sensor 
magnetics [Ref. 3-36], 

5 Sensor characteristics relating the temperature and power sup- 
ply coefficients of offset error and nonlinearity can be de- 
rived empirically. 
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Fig - 3-9 X AXIS NONORTHOGONALITY 


Using small angle approximations, we can solve for the de- 
sired true value of Hx 


HX 1 a 

Hx - Hy Sin e 

(3-la) 

Hx = 

Hx 1 + Hy Sm s 

(3-lb) 


By measuring e (Chapter V)- and storing the angle as a con- 
stant, the x axis data was then restored using equation 3-lb 
above in subroutine "ORTH" . 

2) Subroutines ROTX and ROTY 

These subroutines compute arithmetic values for Hxh and 
Hyh of equatidn 2- lib using sign magnitude quantities and 
table lookup to determine solutions for the transcendental 
functions. Subroutines "SADD" and "SMPY" are nested and used 
to perform double precision add and multiply as required. 
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3) Subroutine HVEC 

Following computation of the horizontal X and Y axis mag- 
netic vector, the subroutine "MAIN" calls subroutine "HVEC" to. 
compute the square of the horizontal vector. Vectors Hx and 
Hy are squared by calling subroutine "SQU" then added, yield- 
ing H (HORIZONTAL) 2 . 

4) Subroutine WTCH 

To compute heading, equation 2-1 (or a similar form) must 
be solved using the horizontal magnetic field vector and either 
the x or y axis horizontal field component. Although the 
square root operation implied in equation 2-1 could be imple- 
mented using a numerical technique [Ref. 3-37, 3-38], the 
computation time is decreased by using a table lookup method. 
Subroutine "WICH" (Fig. 3-10) compares the absolute magnitude 
of the two horizontal field vectors Hx and Hy to determine the 
relative heading of the aircraft 6 with respect to the north- 
south’ and east-west axes (Fig. 3-11) . 



Fig. 3-10. SUBROUTINE "WICH" 


6 If |Hx|<|Hy|, then an equation similar in form to 2-la must 
be used. 
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N 



Fig. 3-11. MAGNITUDES OF Hx AND Hy RELATED 
AIRCRAFT HEADING 

5) Subroutines COSY and SINY (Fig. 3-12, 3-13) 

Depending on the relative absolute magnitudes of Hx and 
Hy, either "COSY" or "SINY" is called to compute aircraft 
heading. These subroutines invoke subroutine "DIVI" to form 
the quotient of the axis vector squared and the horizontal 
field vector squared (a double precision operation) . Sub- 
routine "ANGL" is then called to perform an associative table 
lookup operation using successive approximation and inter- 
polation to complete the inverse cos squared operation. The 
double precision binary quantity is then converted to three 
’ digit binary coded decimal format (BCD) prior to computation 
of aircraft heading (subroutine "HDG" ) . 

The subroutine "SINY" of Fig. 3-12 includes a subtraction 
of 'the computed angle from 90 degrees following conversion to 
BCD format. This operation ensures that the angle passed to 







Fig. 3-13. 


SUBROUTINE "SINY 
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the calling subroutine upon exiting either "SINY" or "COSY" is 
an aircraft heading angle relating sensor x to the north-south 
axis . 

6) Subroutine HDG (Fig. 3-14) 

The function of this subroutine is to compute aircraft 

heading having established the angle between the x axis sen- 

* 

sor and the north- south geodetic axis. Determination of the 
actual heading is accomplished by comparing the signs of both 
the x and y axis horizontal vectors prior to computing heading 
(Fig. 3-15) . It should be noted that all of the preceding 
computations leading to horizontal vector data were on sign 
magnitude quantities preserving the correct horizontal vector 
polarities 7 . 

3-5 CONCLUSIONS 


- This chapter has outlined the practical aspects of design- 
ing an instrument to evaluate both the heading algorithms and 
solid state magnetic indicator proposed in previous chapters . 
The chapter outlined a design approach that can be used to 
implement a microprocessor based instrument. In particular, 
the need to consider the total system hardware requirements 
while simultaneously considering the programming requirements 
was identified. Design proceeded by outlining a system block 
diagram (Fig. 3-1) with major subsystems considered. The in- 
strument required a special purpose computer with an analog 
subsystem to sample and digitize five sensor signals. Timing 
and control of the analog subsystem plus digital processing 
of data was controlled by a microprocessor based central pro- 
cessing unit (CPU) . Memory for permanent storage of 


7 It is possible at certain attitudes to require sign revers- 
als when .computing horizontal vectors . 
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instructions and temporary storage of data was implemented 
using memory chips organized on cards with 2048 byte capacity. 
The particular memory chips selected feature pin compatibil- 
ity 8 with both read only and volatile random access versions. 
System inputs consisted of sensor signals from a three axis 
solid state fluxgate magnetometer plus two analog signals sim- 
ulating gyroscope outputs. System outputs consist of visual 
seven segment readout displaying computed heading. In addi- 
tion, an RS-232 teletype interface was provided to facilitate 
system development and experimentanon. 

By identifying the total system m block diagram form at 
the very beginning, the role and requirements of each subsys- 
tem as well as the supporting software were identified. The 
design then evolved on a modular basis with each subsystem and 
its supporting program developed m parallel. In this manner 
pm assignments for input/output ports and critical timing re- 
quirements that involved both hardware and software considera-. 
tion were handled efficiently. By outlining the program re- 
quirements m flow chart form (analagous to the block diagram 
of the hardware subsystem) , subroutines were identified facil- 
itating a modular program development. Where possible, sub- 
routines were shared in a nested manner avoiding replication of 
programming and waste of memory. 

Details of error analysis and calculation of overall sys- 
tem throughput rate were deferred to Chapter IV. It was point- 
ed out however, that errors induced by imprecision of data plus 
truncation and roundoff during processing of the algorithm were 
to be considered early in the design phase. These data were 
required to select the sensors and the analog to digital 


8 Memory integrated circuit (IC) devices of both types can be 
used in the same mechanical sockets with actual chip type 
being used transparent to the remainder of the system. 
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converter as well as to design the supportive software for 
the analog subsystem. In addition, the data precision require- 
ments were necessary prior to programming, the algorithm 9 . 

By incorporating a microprocessor as the main CPU element, 
considerable sophistication in both control and computing per- 
formance was achieved. The overall system was designed rel- 
atively quickly , provided a convenient laboratory instrument 
for evaluation of the proposed algorithms and featured inhere 
ent flexibility. 


9 Some of the subroutines required double precision manipulations 
to maintain overall system accuracy. 
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CHAPTER IV 

HEADING INSTRUMENT ERROR ANALYSIS 
4-1 INTRODUCTION 

The heading instrument designed to evaluate the heading 
and solid state remote magnetic indicator algorithms is prone 
to error from many sources. These errors will accumulate and 
degrade the accuracy of aircraft heading or yaw angle computa- 
tions. This chapter addresses the various error sources to 
determine -their relative magnitudes and effects on the overall 
computation. 

Prior to beginning the hardware design of the micropro- 
cessor based instrument many of these potential error sources 
were considered. Their effects were considered in. establish- 
ing parameters such as word lengths, A/D converter precision, 
computation speeds, sampling rates, magnetometer sensor accur- 
acies', system noise tolerance, etc. As the design of the 
microprocessor based system evolved, the error analysis re- 
fined. Ultimately, important limitations in instrument: design 
and operation were identified by combined error analysis and 
empirical data. By carefully analyzing the source and extent 
of the limiting parameters (such as sensor offset and non- 
orthogonality) , the magnitude of errors unique to this lab- 
oratory sensor array were identified. Specialized software 
was then added (with empiracally derived constants) to correct 
for the otherwise limiting sensor irregularities improving the 
total system performance. 

In this manner, it is apparent that error analysis is an 
integral part of instrument design. Not only are important para- 
meters identified early in the design cycle (prior to system 
block diagram development) , but shortcomings in conventional 
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sensors can be improved by judicial application of error cor- 
recting algorithms. In this case, data constants were deter- 
mined after the final instrument became operational. The sen- 
sor peculiarities were analyzed empirically using the instru- 
ment itself. 

The chapter begins by first identifying and carefully 
analyzing potential error sources m the sensors. This analy- 
sis is followed by a similar consideration of errors origin- 
ating in the analog subsystem. Processing errors that origin- 
ate due to the finite word length and precision of the micro- 
processor along with che effects of simplifications made to 
the algorithms are finally analyzed. The chapter then con- 
cludes with a summary of measurement errors , a sample error 
analysis, a comparison of predicted to measured error and a 
summary. 

4-2 SENSOR ERRORS 


The heading computation algorithm employing the remote 
magnetic indicator (Chapter II) is prone to error proportional 
to both fluxgate magnetometer’ sensor and gyroscope measurement 
errors. Errors inherent in the fluxgate magnetometer are sum- 
marized on the data sheet [Ref. 3-34] . Since the experiment- 
ation employed simulated gyroscope sensors with voltage levels 
accurately represented, the analysis of sensor errors will 
assume ideal gyroscope sensors to predict experimental data. 

A) Sensor Offset Error 

Magnetometer sensors exhibit error caused by both elec- 
tronic and magnetic phenomena. Errors in the Develco sensors 
were outlined by Workentine [Ref. 4-1] . These offset errors 
are induced in the Develco sensors by both electronic offset 
voltages and currents in the respective sensor electronics and 
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by residual magnetic fields- in the magnetic mass of the sensor 
assemblies. Although the physical and electronic design at- 
tempts to reduce offset error, a finite non-zero output can 
exist when a- zero input is applied. 

Offset error for each sensor in the Develco model 9200C 
three axis magnetometer assembly is specified [Ref. 3-341 as 
"Zero Field Bias +2.5 Volts +1.0%". This offset translates 
into a worst case maximum error voltage of 

EOFFSET = + (2 . 5V x 0.01) = +25mV 

Since the offset error is sensor dependent, correction cannot 
be made at a single physical point (as for analog subsystem 
offsets described in Section 4-3) . Corrections can however be 
made to the measured data by simply adding, or subtracting a 
constant equal to the offset magnitude following each data 
measurement 1 . 

'Offset values for each sensor used in the experiment were 
obtained by rotating the sensor into alignment with earth’s 
magnetic field vector to measure both positive and negative 
maximum values. The difference in magnetic measurement (assum- 
ing negligible analog subsystem error) is related to system 
offset error composed of sensor electronic and sensor plus test 
fixture induced magnetic offset error. The actual offset error 
can be calculated using these two measurements 


Offset corrections were made in the sample subroutine "SAMP" 
illustrated in Fig. 3-6a. 
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| Emax | = Ef + Eo 
| Emin | = Ef - Eo 

| Emax | - [Emm | = (Ef + Eo) - (Ef - Eo) = 2 Eo 
Eo = (1/2) {[Emax | - |,Emin|) 


where 


Emax = The maximum positive voltage recorded when the 
sensor aligns with earth's field vector. 

Emin = The maximum negative voltage recorded t when the 
sensor aligns 180° with earth's field vector. 

Ef = The magnitude of earth's magnetic vector repre- 
sented in volts. 

Eo = The sensor offset voltage due to both electronic 
and magnetic phenomena 

Data 'recorded during x, y and z axis offset measurements as 
described above are recorded m Table 4-1. Since the offset 
error is a function of sensor magnetic permeability, the actual 
offset value will vary with time depending on induced magnetic 
fields 2 . 

Final offset correction values were determined by rotat- 
ing two sensors in the horizontal plane around the third ver- 
tical axis and measuring offsets in two sensors at a time. 
Recorded data for each sensor was previously corrected for 
orthogonality error by the sample subroutine "SAMP" (discussion 


2 For example, magnetized screwdrivers or other tools used near 
the sensor will alter the residual magnetic field. 
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SENSOR 

AXIS 

. DATA RECORDED 
| Emax | 

(HEXADECIMAL) 
[ Emin | 

OFFSET 
(HEXADECIMAL 
PLUS SIGN) 

X 

628 

E78 

+ 40 

Y 

640 

E5D 

-15 

Z 

637 

E68 

-25 


Table 4-1 OFFSET DATA DERIVED BY MEASURING 
EARTH'S FIELD 
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of this correction follows in Section 4— 2B) . Data recorded in 
this manner appears in Tables 4-2 and 4-3. Final correction . 
terms for correcting sensor offset error were calculated using 
these data. Offset terms to be added or subtracted from res- 
pective data channels are tabulated in Table 4-4. 

By correcting system offset errors in this manner, the 
effective error contribution can be reduced appreciably (see 
final data discussion Chapter V) . For a flight instrument, 
sensor offset characteristics as a function of temperature 
variation and supply voltage can be derived empirically and 
appropriate offset corrections made by computing the value of 
the correction term variable . Magnetically induced offsets can 
be reduced by degaussing the sensor assembly periodically. 

B) Axis Alignment: Errors 

The error specification of [Ref. 3 — 34'] indicates' that the 
maximum axis alignment error is +1 degree relative to base ref- 
erenced coordinates. This error results in sensor directional 
uncertainty as illustrated m Fig. 4-1. Each sensor is located 
within a right circular cone with axis along the true sensor 
axis and vertex at the common sensor origin. Although this 
alignment uncertainty contributes no error m determining the 
total magnetic vector 

_ 2 2 _ 2 t 

H = (Hx + Hy + Hz ) , 

there is considerable uncertainty in attempting to resolve the 
true magnetic field component along any axis of the reference 
coordinate system. This alignment uncertainty of magnetic 
sensors limits system performance of conventional field dir- 
ection measuring apparatus [Ref. 4-1] . 



Protractor 

Heading 

Measurement 

(Degrees) 

Data 

Measured 
Hx Hy 
(Units) 

Protractor 

Heading 

Measurement 

(Degrees) 

Data 

Measured 
Hx Hy 
(Units) 

Error Due To 
Offset 
X Y 

(Units) 

0 

8 

-7S9 

180 

83 

720 

75 

-39 

345 

199 

-743 

165 

-112 

704 

87 

-39 

330 

383 

-677 

150 

-292 

638 

91 

-39 

315 

540 

-569 

135 

-448 

527 

92 

-42 

300 

666 

-421 

120 

-573 

381 

93 

-40 

285 

7 4'9 

-249 

105 

-658 

206 

91 

-43 

270 

783 

- 62 

90 

-696 

14 

87 

-48 

255 

770 

132 

75 

-681 

-179 

89 

-47 

240 

706 

314 

60 

-6 16 

-362 

90 

-48 

225 

597 

473 

45 

-504 

-521 

93 

-48 

210 

449 

601 

30 

-359 

-645 

9-0 

-44 

195 

273 

687 







(DATA UNAVAILABLE DUE TO TEXT FIXTURE 
165 LIMITATION) 


TOTAL OFFSETS 


978 477 


AVERAGE OFFSETS 


83.9 43.4 


Table 4-2 X AND Y AXIS ERROR MEASURED BY ROTATING 
X, Y AROUND Z IN THE HORIZONTAL PLANE 
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Protractor 

Heading 

Measurement 

(Degrees) 

Hz Data 
Measured 
(Units) 

Protractor 

Heading 

Measurement 

(Degrees) 

Hz Data 
Measured 
(Units) 

Of f set 
Error 
(Units) 

0.5 

0 

180.5 

-56 

-56 

315.5 

--527 

135.5 

469 

-58 

270.5 

-763 

90.5 

707 

-56 

225.5 

-571 

45.5 

513 

-58 


TOTAL OFFSET -228 

AVERAGE OFFSET -57 


Table 4-3 Z AXIS OFFSET ERROR MEASURED BY ROTATING 
THE S AXIS AROUND THE VERTICAL X AXIS 


Sensor Total Average Required Amount of Correction 

Axis Offset (Units) Correction Decimal Binary Hex 


X 

88.9 

Subtraction 

45 

00101101 

02D0 

Y 

43.4 

Addition 

22 

00010110 

0160 

Z 

47.0 

Addition 

29 

00011101 

01D0 


Table 4-4 OFFSET CORRECTION VALUES 
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Fig. 4-1 SENSOR ALIGNMENT UNCERTAINTY 

Although this error source can be reduced by physically 
aligning the sensors more accurately during assembly, cost of 
the sensors increases. Ultimately, directionality of the mag- 
netic sensors becomes a function of the physical sensor itself 
and more accurate sensors are required as pointed out by Gise 
[Ref. 4-2] . A heading system that tolerates sensor misalign- 
ment is therefore a very desirable alternative to requiring 
precise alignment or more elaborate sensors. 

During assembly of the Develco fluxgate magnetometer sen- 
sor array, sensor misalignment is determined by using earth's 
magnetic field and a precision mechanical rotation assembly. 

A sensor (assume the X axis) is aligned with earth's magnetic 
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vector by positioning the sensor to maximize electrical output 3 . 
One of the other sensors (assume the y axis-} is aligned with the 
rotation axis of the precision calibration assembly (Fig. 4-2) 
and perpendicular to the first by rotating the sensor array 
around the second sensor axis (y axis m this case) and adjust-, 
ing its relative position until a null output is achieved at all 
rotation angles. Mechanical orthogonality or the sensors is then 
limited only by the mechanical imprecision of the calibration 
device (orthogonality within tO.Ol degrees can be eaaily achieved 
i-n the calibration tool) and by the directional characteristics 
of the physical sensors. 


X Axis 



Fig. 4-2 MECHANICAL ORIENTATION OF THE MAGNETOMETER 
SENSORS DURING CALIBRATION 

In addition to functioning as an alignment apparatus , the 
calibration device described above provides a convenient means 


3 By maximizing or nulling a measurement, the mechanical posi- 
tioning is a function of only the field and the resolution of 
the voltage measuring device obviating errors due to physical 
position measurement. 



64 


to characterize sensor assemblies after final assembly adjust- 
ments are made- Any misalignment of the second sensor relat- 
ive to the first results in a coning of the second sensor 
around the rotation axis 4 with a sinusoidal output voltage 
that is a function of total earth's magnetic field and axis 
alignment error. The peak to peak voltage resulting from sen- 
sor coning is recorded during the' final alignment test and 
made available to sensor purchasers. Coning voltages developed 
for the sensor assembly used with this experiment were obtain- 
ed from Develco [Ref. 4-3] and are recorded in Table 4-5. Sen- 
sor misalignment for each axis can be derived using additional 
data provided by Develco along with additional empirical data 
derived by experimentation. 

The total ambient magnetic field at the Develco laboratory 
is measured using the three sensors (applying equation 4-1) 
and is supplied as digital data. In our case, the total field 
measured was 15.73 units or 


• 1573 units 

2048 units 
Pull Scale 


x 60,000 
(F.S.) 


gamma F.S. 


46,084 gamma (y) 


Sensitivity of the sensor = - 1 = 42 piVolts/y 

Considering the X axis sensor, coning resulted in a signal of 
38 mV peak to peak (or 19mV peak) . Misalignment of the X axis 
sensor from the Y-Z plane can then be calculated as 


4 Assume that the first axis is initially adjusted for maximum 
output to align it with earth's field and the rotation axis 
is perpendicular to the field. 
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Peak Signal = 19 mV or 456 gamma angular misalignment 

e.. = Sin" 1 4‘56 
46,084 

e = 0.57 degrees 

A 

Similarly, the Y and Z axis have misalignment errors of 
ey - 0 and ez = 0.76 degrees with respect to the X-Z and X-Y 
planes respectively (sensor orthogonality errors are tabulated 
in Table 4-5) . 

Having established that sensor orthogonality errors exist, 
the remaining task is to identify the direction that the sensor 
axis points relative to the other two sensor axes. Since the 

Y axis has relatively little orthogonality error’, it will be 
assumed to be perpendicular to the X-Z plane. In addition, 
since the Hz data enters into the algorithm in a second order 
manner relative to the Hx and Hy measured data, correction and 
characterization of the Hx sensor was considered to be of prim- 
ary concern. Orientation of the X axis sensor relative to the 

Y and Z axes was determined empirically. 

Angular position of the X axis sensor can be described 
using the error angles £xy and exz as delineated in Fig. 4-3. 
Characterization of sensor orthogonality error m terms of 
these two angles would enable algorithmic corrections of mea- 
sured data. 



Fig. 4-3 X AXIS SENSOR ORIENTATION 
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1) Empirical Determination of exz 

The .angle exz (angle between the x axis sensor and the z. 
axis of the geodetic coordinate system) was determined in sev- 
eral steps using the test apparatus described in Chapter V. 

i) The x and y sensors were oriented m the horizontal 
plane with the z axis sensor vertical downward. 

ii) The x and y sensors were rotated around the z axis 

with magnetic data measurements (corrected for sensor 
offset error as described in Section 4-2A) recorded 
in Table 4-6 for incremental rotation angles. 

in) The total horizontal field at each angular position 
was calculated 

Hht = <Hx 2 + Hy 2 )*. 

iv) Average horizontal field Hav was computed by averag- 
ing the results of iii) above. 

v) The horizontal field deviation Hd was computed for 
each angular position; tabulated in Table 4-6 and 
plotted on Fig. 4-4. 

Hd = (Hav - Hht) 

The horizontal field deviation or error (as shown on Fig. 
4-4) was now examined. An angular error exz should cause the 
horizontal field error curve to peak at 90 and 180 degrees. 
Since this obviously was not the case, it was concluded that 
major error m x axis orthogonality was due to the component 
cxy . 
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Total 

Computed 


Physical* 

Heading 

Displayed** 

Heading 

Measured Data 
(Units) 

Horizontal 
Field (Hht) 

Hd 

(Hav-Hht) 

(Degrees) 

(Degrees) 

Hx 

Hy 

(Units) 

(Units) 

355 

90 

12 • 

-727 

727 

-3 

335 

70 

258 

-682 

729 

-1 

315 

50 

477 

.-555 

732 

2 

295 

30 

637 

-362 

733 

3 

275 

10 

720 

-127 

731 

1 

255 

350 

718 

12 7 

729 

-1 

235 

330 

628 

363 

725 

-5 

215 

310 

454 

569 

734 

4 

195 

290 

242 

687 

728 

-2 

175 

270 

-10 

731 

731 

1 

155 

250 

-258 

686 

733 

3 

135 

230 

-476 

557 

733 

3 

115 

210 

-635 

362 

731 

1 

95 

190 

-721 

123 

731 

1 

' 75 

170 

-718 

-133 

730 

0 

55 

150 

-628 

-371 

729 

-1 

35 

130 

-460 

-565 

729 

-1- 



Total Hht = 12415 




Average 

(Hav) = 

730 


*Measured 

using a protractor 

on the 

test apparatus 



**Computed and displayed digitally by the instrument. 


Table 4-6 MEASUREMENT OP HORIZONTAL FIELD 




Fig. 4-4 


DEVIATIONS OF THE HORIZONTAL FIELD MEASUREMENT 
FROM THE MEAN 


CTk 

LO 





71 


2) Empirical Determination of sxy 

The angle exy representing x axis sensor misalignment 
relative to axis y was measured as follows : 

Steps i) and ii) above were repeated with the exception 
that the calculated values for Hx and Hy (Hxc and Hyc respect- 
ively) were recorded with measured Hx and Hy data (Hxm and Hym 
respectively) in Table 4-7. The calculated values were obtain- 
ed by assuming that the angle exz as determined above was neg- 
ligible and that the y axis sensor was perpendicular to the x-z 
plane. With these assumptions, we note that at the heading of 
zero degrees (extrapolated between display of 10 and 350 degrees 
of Table 4-7 and Fig. 4-5) , there is no error in yaw due to 
either Hx or Hy. By physically rotating the sensors in fixed 
intervals from yaw = 0 degrees and noting that the horizontal 
field Hh = 730 units, we can then compute expected Hx and Hy 
data at respective yaw orientations. 

Physical orientation of the x axis sensor is easily deter- 
mined by considering orientation at the maximum error excurs- 
ions. These observations are illustrated m Fig. 4-6. We note 
that the only possible orientation of the x axis sensor satis^ 
fying the data m Fig. 4-5 is that of Fig. 4-6. 




Fig. 4-6 (a) SENSORS ORIENTED AT YAW = +90 degrees 
(b) SENSORS ORIENTED AT YAW = +270 degrees 



* 

** 


Physical* 

Heading 

(Degrees) 

Displayed** 

Heading 

(Degrees) 

Measured Data 
(Units) 

Hxm Hym 

Computed Data 
(Units 
Hxc Hyc 

Deviation 

(Hxm-Hxc) 

(Units) 

(Hym-Hy 

355 

90 

12 

-727 

0 

-730 

+12 

+ 3 

335 

70 

258 

-682 

250 

— 686 

+ 8 

+ 4 

315 

50 

477 

-555 

469 

-559 

+ 8 

+ 4 

295 

30 

637 

-362 

632 

-365 

+ 5 

+ 3 

275 

10 

720 

-127 

719 

-127 

+ 1 

0 

255 

350 

718 

127 

719 

+127 

- 1 

0 

235 

330 

628 

363 

632 

365 

- 4 

- 2 

215 

310 

464 

569 

469 

559 

- 5 

lp 

195 

290 

242 

687 

250 

686 

- 8 

1 

175 

270 

-10 

731 

0 

730 

-10 

1 

155 

250 

-258 

686 

-250 

686 

- 8 

0 

135 

230 

-476 

557 

-469 

559 

- 7 

- 2 

115 

210 

-635 

362 

-632 

365 

- 3 

- 3 

95 

190 

-721 

123 

-719 

+127 

- 2 

- 4 

75 

170 

-718 

-133 

-719 

-127 

1 

- 6 

55 

150 

-628 

-371 

-632 

-365 

4 

- 6 

35 

130 

-460 

-565 

-489 

-559 

9 

- 6 


Measured using a protractor on the test apparatus. 
Computed and displayed digitally by the instrument. 


Table 4-7 MEASURED AND COMPUTED Hx AND Hy DATA 
IN THE HORIZONTAL PLANE 



73 


Magnitude' of the angle exy can be“ computed as follows 
using data from Fig. 4-5 

Max. delta from Fig. 4-5 = 10 units 
Average horizontal field = 730 units 

exy max = Sin ^ yyg- 

= 0.79 degrees 

We note that the angle of 0.79 degrees is approximately 
the same as determined by Develco during manufacture of the 
sensors {Table 4-5) . The added error is due to test set in- 
accuracy. 

C) Fluxgate Sensor Noise Induced Error 

The analog output from the fluxgate sensors can exhibit an 
error due to signal uncertainty resulting from noise. Although 
the data sheet [Ref. 3-34] indicates that 5mV peak to peak of 
ripple can exist on the output, the frequency content centers 
m the 550 kHz range (driver frequency of the fluxgate magneto- 
meter) and no appreciable ripple 5 exists below 60 Hz (especial- 
ly when the sensor output is filtered prior to data sampling) . 
The noise specification of less than 1 gamma peak to peak in 
the 1 Hz bandwidth region is also negligible. In summary, no 
appreciable error due to noise on the magnetometer signal lines 
is evident. 

D) Magnetometer Gain Error 

The magnetometer is specified to have gain (sensitivity) 
of 2,. 5 Volts/600 milligauss, +1% which translates into a maxi- 
mum signal uncertainty of 

5 Verbally confirmed by Workentine of Develco [Ref. 4-1] . 
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+ ( 2 . 5V x 0.01) = +25 mV. 

This represents a sensor transfer function of 4.16 Volts/gauss 
or 0.24 gauss per volt. The uncertainty then can be expressed 
as 


+(0.24 gauss x 0.01) = +2.4 mxlligauss 

- = +(2.4 x IQ 2 ) gamma 

Since this error is not corrected in the laboratory instrument 
it. will be considered in total m the final error analysis. It 
is worth noting however, that should the magnetometer gain un- 
certainty be characterized, gain corrections for each sensor 
could be made during computation by the computer. In addition 
the error term is proportional to actual signal level applied. 

E) Magnetometer Linearity Error 

D.C. linearity of the magnetometer is specified to be 
+0.5% of signal level. This uncertainty at full scale can be 
expressed as +(2.5 Volts x 0.005) = +12.5mV. Alternately, lin- 
earity error can cause a signal uncertainty of +1.2 milligauss 
or +(1.2 x 10 2 ) gamma. Linearity error is also nor corrected 
during computation and is considered m the final error analysi 
By simply characterizing and correcting the -linearity character 
istics of each sensor, considerable improvement in system accur 
acy could be achieved. 

4-3 ANALOG SUBSYSTEM ERROR ANALYSIS 

The analog subsystem of the instrument is outlined in 
block diagram form in Fig. 4-7. This subsystem accepts . analog 
signals from magnetometer and gyroscope transducers , performs a 
time division multiplexing between the signals and digitizes 



+5 


+15 


-15 



Fig . 4-7 THE ANALOG SUBSYSTEM 


ui 
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the respective signals prior to subsequent processing by the 
computer. During this data acquisition and conversion process, 
errors are introduced into each of the signals. This section 
addresses the potential error sources and computes the respect- 
ive error contributions to be expected during operation of the 
instrument. 

Although the multiplexer and sample and hold blocks of Fig. 

.4-7 could be eliminated (eliminating possible error sources) by 
digitizing each signal with a unique analog to digital converter, 
it can be shown that such a system would be expensive and dif- 
ficult to implement. The analog to digital converter (A/D) 
quantizes an analog signal in a finite amount of time. Speed of 
conversion is predicted in a finite amount of time by both the 
resolution of the converter and 'the frequency of the signal to 
be converted. Time required to perform a conversion is general- 
ly called the "aperature time". 



Fig. 4-8 APERATURE TIME AND AMPLITUDE UNCERTAINTY 

As illustrated in Fig. 4-8, aperature time and amplitude 
uncertainty are related by the time rate of change of the ana- 
log signal. For the particular case of a sinusoidal signal to 
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be converted, the maximum rate of change occurs at the zero cross- 
ing of the waveform and the amplitude change, is : 


” A 

AV = dt t = 0 x ta {4-la) 

AV = V w ta (4-lb) 

giving ~ = w ta = 2 irf ta. (4-2) 

From this result we can determine the aperature rime required to 
digitize a 30 Hz signal to 12 bits resolution (a resolution of 1 
part in 4096 or 0.0244%). 


ta 


V „ 1 .000244 

V X 2 f 6.28 x 30 


1.3 x 10 


This result indicates that to remain within 1 bit of resolution 
(0.0244%) we require an aperature time of 1.3 microseconds to pro- 
cess analog signals varying at a rate of 30 Hertz. It can be seen 
that the system would require fast A/D converters plus extremely 
fast computational capability to accommodate this configuration of 
sensors and analog subsystem. By using multiplexing and sample 
and hold circuitry we can however reduce the number of A/D con- 
verters required to one and alleviate the aperature and processing 
requirements imposed above. 

The operation of sampling to be used by the instrument is 
illustrated in Fig. 4-9 which shows an analog signal and a train 
of sampling pulses. The pulses are provided by the central pro- 
cessing unit. A switch connects the analog signal for a very 
short period of time to the hold circuitry charging a capacitor 
and storing the sampled voltage until the next sample is required. 
This type of sampler is called sample and hold. 
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SIGNAL 


SAMPLING PULSES 


SAMPLED SIGNAL 


SAMPLED AND 
HELD SIGNAL 


Fig. 4-9 SIGNAL SAMPLING PROCESS 
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A) Sampling Rate Errors 

The process of uniformly sampling a function of contin- 
uous time can yield a significant source of error if the sampl- 
ing period T is .selected too large [Ref. 4-4, 4-5]. This error 
can be illustrated by considering an analog signal xa(t) that 
has the Fourier representation [Ref. 4-6] 


xa(t) = 2 ~ a(3^)e^ t d$2 

(4-3a) 

Xa{]^) =■ r_“xa (t) e"^ t dt 

( 4-3b) 


The sequence x(n) with values x(n) = xa(nT) is said to be de- 
rived from xa(t) by periodic sampling and T is the sampling 
period. The reciprocal of T is called the sampling frequency 
or' sampling rate. In order to determine the sense in which 
x(n) represents the original signal xa(t) , it is convenient to 
relate Xa(jSl), the continuous-time Fourier transform of xa(t), 
to X(e J ), the discrete-time Fourier transform of the sequence 
x(n). From (4-3a) we note that 

x(n) = xa(nt) = ~ a(Dft) e^^dft (4-4) 

From the discrete-time Fourier transform we also obtain the 
representation [Ref. 4-4] 


x (n) = ^ 


(4-5) 


To relate the equations (4-4) and (4-5) we can express (4-4) as 
a sum of integrals over intervals of length 2it/T, as in 


_1 

2ti 


09 f( 2r+l) tt/T 

X I Xa(jfi) 

r=-«J (2r-l) t/T 


e jJmT dfi 


x (n) 


(4-6) 





80 


“Each term in the sum can be reduced to an integral over the 
■•ange -ir/T to +n/T, by a change of variables to obtain 


x(n) = tT ^ 


tt/T 


Xa 


r=-« J -tt/T 


“ r tt/T 

(nl-jiD / Xatja+j2^)e :)0n V 2,rrn dQ 
r»— J -ir/T 1 


j (fl+^) 


eJ (Q + ^£)nT dn 


( 4-7a) 


x 


(4-7b) 


I f we now change the order of integration and summation and 
~i 2 tt rn 

note that e =1 for all integer values of r and n, we ob- 


tain 


x(n) = 


2it 


ir/T 

-ir/T 


S Xa(jJ2'+j^0 

r=— * 


e J dfi 


(4-8) 


|y substituting Q = u)/T we get 


x(n) = ji 


-TT 


b ^ xa{i^4j^) 


T — 'T 
r=-= 


nun, 
e J dw 


(4-9) 


— 1 — 

Rhich is identical ,xn form to equation (4-5). We can therefore 
jake the identification (equating like terms of (4-5) and (4-9) 

oo 

X (e-^) = i E xa(^4 3 ^) (4-10) 


r=-«° 


J© can also express (4-10) in terms of the analog frequency var- 
iable fi (where tt - oj/T) as 


X (e^ uT ) = ~ £ Xa ( jQ+j-—^-) 


(4-11) 


r =- c 
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The last two equations clearly reveal the relationship 
between the continuous -time Fourier transform and the Fourier 

i 

transform of a sequence derived by sampling. For example, if 
Xa (jfi). is as depicted in Fig. 4-l'0a then XCe^) will be as 
shown in Fig. 4-10b when the sampling period T is too long and 
as shown in Fig. 4-1QC if T is short enough. 

. ooT % 

From Fig. 4-10c it is obvious that if — <ir, i.e., we 

sample at a rate at least twice the highest frequency of Xa(jft) , 
then X(e-* w ) is identical to Xa(u/T) in the interval -ir<w<_n and 
can be recovered from the samples xa(nT) by an appropriate 
interpolation formula. 

For the remote magnetic indicator instrument designed in* 
previous chapters, rhe analog signals are filtered with a low 
pass sec-cion reducing frequency content above 30 Hz. The sampl- 
ing rate must therefore exceed 60 Hz (T<16.67 m.s.) to enable 
accurate dynamic operation of the system. 

Laboratory measurements of sampling rates on the function- 
al microprocessor based instrument revealed that the analog 
subsystem operated at a sampling rate of 62.5 Hz (16 m.s.) in- 
dicating that the algorithm execution rate supported a system 
bandwidth of 31.25 Hz. If frequency content of the analog 
signals is less than 31.25 Hz there is no error due to sampl- 
ing . 

B) Analog Multiplexer Induced Error 

The analog multiplexer of Fig. 4-7 selectively connects 
one analog transducer output at a time to the input of the 
sample and hold subsystem. The Datel Systems, Inc., multiplexer 
(Ref. 4-7] selected for the remote magnetic indicator experiment 
features eight MOS-FET switches with associated driver circuits , 



Xa ( j R ' 

\ 


(a) Fourier transform as a 

continuous -time signal 

- - 

X(e^ W ) 




Fourier transform of the discrete- time signal obtained by 
periodic sampling (T is too large) 


(c) Same as (b) except T is short enough 


Fig. 4-10 FOURIER TRANSFORMS OF CONTINUOUS AND 
DISCRETE-TIME SIGNALS 
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FET pull-up to reduce propogation delays and all of the nec- 
essary decoding logic to enable random channel addressing with 
a four bit parallel binary input. 

Several important parameters are used to characterize ana- 
log multiplexers and can contribute error. 

1) Transfer Accuracy 

Transfer accuracy is a function of the source impedance, 

switch resistance, load impedance (if the multiplexer is not 

buffered) and the signal frequency. It expresses the input to 

output error as a percentage of the input. In our case the 

* 

system configuration predicates a maximum error due to transfer 
accuracy of (+0.01%) yielding an error term of 

+0.0001 x 2.5 Volts = +25 mV 


2) Settling Time 

This parameter defines the time elapsed from the applica- 
tion of a full scale step input to the time when the output has 
entered and remained within a specified error band around its 
final value. In our case rhe selected multiplexer has a maxi- 
mum settling time of 1 microsecond to +0.01% full scale (F.S.) 
Since the control system selecting channels is implemented 
using a microprocessor, the minimum time between analog sub- 
system commands will always be greater than 3.0 microsecond 6 . 

The multiplexer will therefore always have settled to the final 
value before the sample and hold circuit (following this sub- 
system) can be activated with no error due to the settling time 
parameter . 


6 One machine cycle, time for the 2650 microprocessor with 1 mHz 
clock frequency. 
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3) Throughput Rate 

The highest rate at which the multiplexer can switch from 
channel to channel at its specified accuracy is in this case 
500 kHz. Since this rate is more than four orders of magnitude 
greater than the operational rate of the subsystem there is no 
error due to' throughput rate limitations. 

4) Input Leakage Current 

The amount of signal coupled to the output as a percentage 
of input sxgnal applied to all OFF channels together can be cal- 
culated by considering the maximum leakage current specified 
from OFF channels to the ON channel. In our case the maximum 
error signal can be calculated 

Error = [4 (8 na x 2000 ohms source imped. ] T 

Error = 32 microvolts 

Note that m this case the voltage levels are statistically 
independent allowing an R.S.S. of error sources to calculate 
total error [Ref. 4-8, 4-9] . 

C) Sample and Hold Circuit Induced Errors 

The sample and. hold subsystem consists of a switch and 
capacitor arrangement as shown in Fig. 4-11. The Datel Systems, 
Inc., model SHM-IC-1 integrated circuit sample and hold device 
[Ref. 4-10] features a self-contained high gain differential 
input amplifier, a digitally controlled electronic switch and 
a high input impedance buffer amplifier. The external compon- 
ents used with the sample and hold- circuit in’ the solid state 
remote .magnetic indicator instrument consisted of the O.OOlyf 
holding capacitor and a 100K offset tnmpot. By connecting 



he output back to the negative input of the input amplifier 
Fig. 4-11 ) , the sample and hold subsystem operated in a unit 
ain, noninverting mode. When the switch is closed, the unit 
s in the sampling or tracking mode (Digital Control = 0 Volt 
nd will follow a changing input signal. 


— *r« •— — 1 


Input 


Outpui 


11 12 


100K 


Guard Ring 
: 0 . 001 


w\ 
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When the switch opens the unit is m the hold mode and 
retains a voltage on the capacitor for some period of time 
depending on capacitor and switch leakage. Sample and hold 
devices are characterized by a number of important parameters 
that must be considered in the design of a data acquisition 
subsystem. 

1) Acquisition Time 

The time lapse between the time that the sample command 
is given to the point where the output enters and remains with- 
in a specified error band around the input value is specified 
to be less than 4 microseconds time to transit from 0 to 0.1% 
of 10 Volts with C = 0.001 yf [Ref. 4-10] . This implies that 
the control signals emanating from the central processor should 
allow at least 4 ys acquisition time prior to entering the hold 
mode. We note that the sample and hold subroutine (Appendix 
B) executes the instruction 

IORI , R3 H 1 80 ’ READY TO HOLD, 

a two machine cycle instruction prior to sending the hold con- 
trol signal. This instruction delays control signal transmis- 
sion by (2 x 3 ys) = 6 ys allowing the sample and hold circuit 
ample time to settle with no appreciable error due to the 
acquisition time parameter. 

2) Hold Mode Voltage Droop 

The maximum change in output voltage as a function of time 
is specified to be 50 mv/sec maximum using a 0.001 yf poly- 
styrene capacitor. Since the maximum total accumulated time to 
completion of the analog to digital conversion can be calculated 
as 
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5 Instructions (11 machine cycles)* = 33 ys 

1 Analog to Digital Conversion = 20 ys 

3 Instructions if Conversion not 
synchronized with instructions 
(7 machine cycles) = 2_1 ys 

7 4 ys 

_ 5 

we can then compute droop error to be 50 mv/sec x (74 x 10 ) 

sec = 3.73 mv. 

3) Aperature Delay 

The maximum time lapse between the time of hold signal 
receipt to opening of the switch is specified co be 50 nsec, 
an insignificant length of time m the instrument. There is 
therefore no error due to aperature delay. 

4) Offset Error 

Although the maximum offset error is specified to be 20 mv 
maximum [Ref. 4-10] , the error was eliminated using the- 100K 
trimpot offset adjustment. There was no appreciable offset 
error contribution due to the sample and hold circuit. 

5) Gain Error 

The gam error of a sample and hold circuit is apparent 
during the sample mode when the transfer function of the total 
amplifier deviates from the ideal unity slope condition (Fig. 
4-12) . In the noninverting unity gam mode, the specified gam 
error is +0.05% maximum* yielding a signal error of 

(+0.0005) x 5.0 V = 250 mV 

This error can, however, be eliminated with the gain 
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Fig. 4-12 GAIN, OFFSET AND LINEARITY ERRORS 

adjustment available at the analog to digital converter. There 
will ’therefore be no appreciable net gam error due to rhe an- 
alog subsystem. 

6) Nonlinearity Error 

Nonlinearity error is apparent m the sample and hold 
circuit if the transfer function departs from a linear curve 
(Fig. 4-12). In the nonmverring unity gam mode with a 
0.001 yf holding capacitor the maximum nonlinearity is 0.01% 
resulting m a worst case signal uncertainty of (.0001) x 
2.5 Vol*ts = 25 mV. 

7) Hold Mode Feedthrough 

This error appears due to input signal appearing at the 
output when the unit is in the hold mode. Although the 
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feedthrough varies with signal frequency and the expected signal 
frequencies are substantially lower than the upper frequency 
limits of the sample and hold device (30 Hz max. versus several 
kiloHertz) , we consider the worst case feedthrough of 0.01% 

[Ref. 4-10] or 25 mV. 

D) Analog to Digital Converter Induced Errors 

The A/D Converter selected for the solid state magnetic 
indicator instrument (Datel ADC -MAI 2 B IB ) [Ref. 4-11] uses the 
successive approximation technique to achieve excellent lin- 
earity and speed. Important parameters that potentially con- 
tribute errors are addressed below. 

1) Resolution Error 


The smallest analog change that can be distinguished by 
the A/D converter is 


Least Significant Bit (LSB) 


Full Scale 
-,n 


LSB = -=- = 1.22 mV 
2 1 2 


this uncertainty manifests itself as. an error m computing by 
limiting the precision of any calculation. 

2) Linearity Error 

The maximum deviation from a straight line drawn between 
the end points of the converter transfer function are specified 
in [Ref. 4-11] to be + 1/2 LSB (in our case +1.22 mV of analog 
signal) . 
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3) Accuracy Error 

The input to output error of the A/D converter is specif- 
ied in [Ref. 4-11] to be +0.012% F.S. +1/2 LSB or 

+(0.00012) x 5.00V +1.22 mV = +1.82 mV Worst Case 

In reality, the two error terms are unrelated and the 

Rss Error = + [(.00012 x 5V) 2 + (1.22mV) 2 ]* 

RSS Error = +1.36 mV 

4) Offset Error and Gain Error 

Both the offset error and gam error were adjusted to zero 
using the trimming potentiometers (Fig. 4-7) and the calibra- 
tion procedure outlined in Ref. 4-11. A reference signal of 
plus 1/2 LSB (1.22 mV) was applied to -che converter and the 
offset trimming potentiometer adjusted until the output flick- 
ered equally between logic "0" and logic "1". The gain was 
then adjusted by setting the converter input to full scale 
minus 1-1/2 LSB (4.99817 Volts) and the gam trimming potentio- 
meter was adjusted until the output flickered between logic 
"111... 110" and logic "111... 111". The above steps were re- 
peated until no appreciable error in gain or offset was evi- 
dent. 


4-4 PROCESSING ERRORS 

Errors m processing data accrue due to several sources 
including imprecision and truncation. Since the microprocessor 
selected for the instrument is inherently an eight bit device, 
single precision calculations are conducted wiuh eight bits 
and double precision calculations are conducted with a total 
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jof sixteen bits. This section addresses the effects of com- 
iputational precision and truncation in the various subroutines 
and relates these to overall computational accuracy. The var- 
ious subroutines are analyzed in chronological order as they 
appear in the mam program. 

A) Subroutine ”SAMP" 

The sample subroutine (delineated in Fig. 3-6a) selects 
and digitizes analog signals by controlling respective analog 
subsystem modules. During the first portion of this subrout- 
ine, A/D converter data bits are stored m two consecutive 
bytes 7 in the computer memory. The A/D conversion precision 
of 12 bits is thereby preserved. 

The second, third and fourth operations of the sample sub- 
routine convert the unipolar binary format of the data to sign 
magnitude format, adds offset quantities and merely changes the 
[signs of the Hx and Hy data. The operations are conducted m a 
double precision manner and precision of the data remains un- 
altered. 

| Correction of x axis orthogonality error is the final operation 

i 

of the sample subroutine. Equation (3-1) is implemented at 
this .point using a table lookup (for the sin function) , multi- 
plication and addition. The final result can be expressed as 

Hx = Hx' + Hy Sin e 


7 A byte is accepted terminology for an- eight bit data quantity. 
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where the respective quantities have the following forms 


1 1 


Hx J 


- xi + Xy a 2^ 

* . . -i -L 


Hy = 


Sin e = 


i=l 
1 1 

x 2 + X! a 2- 

9=1 J 


£ V k 


and 

X! , x 2 are sign bits 

a^ ^ equal 0 or 1 depending on whether the respective 
term is to exist or not 


We can analyze the effects- of imprecision and truncation 
by noting that the sm function has eight significant binary 
bits resulting in a resolution of 1/256 or 90°/256 = 0.352°. 


The relative error in sin e is computed by Dahlquist 
[Ref. 4-12] as follows 

let a = the approximate value of sm e 
a = the exact value of sm e 

then the relative error in a is 

(a - a) /a if a / 0 

Since data m the sin table has been truncated, maximum 
relative error can be as large as +(1/2 12 ) or +0.02%. 


From the definition of relative error we obtain the 
following relationships between exact, estimate and estimated 
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relative error 

a = a + ar = a Cl + r) 

If a x , and a 2 have relative errors of +0.39% and +0.02%, res- 
pectively, then 

a x a 2 = a x (1 +0.0039) a (1 +0.00024) 

= a x a 2 (1 +0.0039) (1 +0.00024) 


Thus, the relative error in a x a 2 is 

(1 +0.0039) (1 +0.00024) - 1 = 

+ (0.0039) +(0.0039) (0.00024) +(0.00024') 

= +(0.0041) 

Since the maximum value of Sin e to be encountered occurs 
when the orthogonality error (s) is 1 'degree, sin e = 0.017 
maximum. The maximum value for Hy can be 0.6 gauss or 2048 
units. Maximum error due to imprecision in the product HySme 
is then 

Er- Max = (2048 x 0.017) (1 + . 0041) - (2048 x 0.017) 

= 0.1427 units 

Since only the integer portion is retained in the final 
product, insignificant error can be attributed to imprecision 
of the sm e term in this case. Orthogonality error will be 
adequately corrected. 

B) Subroutines ROTX and ROTY 

These subroutines were developed in Chapter III and 
implement the equation of 2-11 required to compute horizontal 
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x and y magnetic field components. Equations to be implemented 
by the respective subroutines are 

Hx = Hxm Cos (pitch) + Hym Sin (pitch) Sm (roll) 

+ Hzm Sin (pitch) Cos (roll) (4-12) 


and 


Hy = Hym Cos (roll) - Hzm Sin (roll) (4-13) 

where Hxm, Hym and Hzm are measured field components made avail- 
able - from the magnetometer via the analog subsystem. 

Since the transcendental functions are implemented using 
cable lookup and are limited in precision to 8 bits, imprecision 
in these variables will dominate in generating error. In part- 
icular, the sin/cos terms will have relative error in the order 
of +1/256 or +0.39% while the measured field , data has relative 
uncertainty of only +1/4096 or +0.02%. Multiplications will 
result m addition of the bounds for the relative error as 
illustrated m section 4-4A above. 

The transcendental terms above are limited in magnitude 
to 1.0 maximum while the field measurements can be 0.60 gauss 
max. In this case the individual product terms of (4-12) and 
(4-13) can have maximum errors of 

Er = (2048) (1 + 0.0041) - (2048) = 8.4 units 

Errors in Hy and Hx (4-12 and 4-13) will be. maximum when 
roll and pitch are at 45 degrees and the fields are equal. In 
this case the error m Hy will be 
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EHy = [{0. 707) (2048) (.1 + 0.0041) - (0. 707) (2048) ] - 
[ (0.707) (2048) (1 - 0.0041) - (0. 707) (2048) ] 

EHy = 4.94 - 5.94 = 11.87 units 

Similarly, maximum error in Hx can be calculated as 

EHx = [(0.707) (2048) (1.0041) - (0 . 707) (2048) ] x3 
EHx = 17.8 units maximum 

It should be noted that these error terms are worst case 
and peak at multiples of 45 degrees m yaw. 

C) Subroutines COSY and SINY 

These two subroutines compute the angle between the x axis 
sensor (when projected onto the horizontal plane) and the 
north-south horizontal vector of earth's magnetic field. The 
first two opeations of these subroutines perform double pre- 
cision multiplication and division. Since the dara variables 
involved are 12 bits in length and the computations performed 
preserving 16 bits, no error is introduced. 

The "ANGL" subroutine called by the above two subroutines 

computes the desired (x axis to horizontal vector) angle by 

completing an associative table look up procedure. The cask 

required is to match a given data quantity either (Hx^/Rh^ or 
2 2 

Hy /Hh ) with the contents of a memory cell. The address of 
this cell is then the required angle. 

•Since the table, is limited in precision to 16 bits there 
are obviously cases where an interpolation is required to 
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ascertain the true address 9 . The function stored in tabular 
form is cos 2 * © where 0 varies from 45 to 90 degrees. Maximum 
error will therefore be induced while attempting to locate 
solutions (angles near 90 degrees if inadequate precision is 
provided. Error in this region due to resolution of tabular 
data can be examined by noting the entries m Table 4-8 


0 

Cos 2 © 

Most Significant 
Binary Bit (2“ x ) 

90 

0 


89 

0.000305 

12 

88 

0.. 001218 

9' 

, 87 

0.00274 

8 


Table 4-8 Cos 2 © AND MOST SIGNIFICANT BINARY DIGITS 

provided to indicate the relative magnitudes of Cos 2 © in the 
region of © = 90 degrees. We observe that the most signifi- 
cant binary digit affected at: 89 degrees is binary decimal 
digit 12 implying that the resolution of Hx 2 /Hh 2 or Hy 2 /Hh 2 
(the argument of Cos 2 ©) must be accurate to at least 1/2 12 or 
0.024%. 


Considering the horizontal field of earth's magnetic 
vector as observed m laboratory experimentation at this lat- 
itude, we note that Hh is 730 units. At a heading of 89 deg- 
rees, Hx = 730 Cos 89 = 12.7 units. The argument would there- 
fore be 

ARG = Hx 2 /Hh 2 a 0.000305 


9 The procedure determines the relative address by linear 

interpolation, then selects the closest address as the 

required angle for the solution. 
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Since, the squaring and division operations are conducted 
in double precision, precision is preserved and the algorithm 
should be able to resolve heading to at least one degree over 
all portions of *the compass . 

D) Errors Due to the Remaining Subroutines 

Since all of the remaining subroutines work with data that 
has been rounded to a precision representing 1 degree or better 
and the computations involve addition or subtraction in double 
precision binary or binary coded decimal (BCD) format, we note 
that there will -be no further appreciable error due "to trunca- 
tion or rounding. 


4-5 MEASUREMENT ERROR SUMMARY 

Errors due to sensors and measurement of their respective 
outputs were discussed in sections 4-2 and 4-3 above. Before 
proceeding with the analysis of errors, the total signal in- 
accuracy due to contribution from the many sources above will 
be summarized in Table 4-10. Total instrument error can then 
be computed by considering the propagation and enhancement of 
these errors during the computation process. 

Since the errors m Table 4-9 are stochastically independ- 
ent, we can compute error for any given signal level by find- 
ing the RSS of respective error sources. In this manner, the 
instrument error can be evaluated by considering all input; 
signals with errors superimposed' to produce an erroneous com- 
putation of heading. 
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PARAMETER 

ERROR 

COMMENT 

Magnetometer 

Offset 

= 0 

Corrected by software 

Orthogonality 

= 0 

Corrected by software 

Noise 

Negligible 


Gain 

+0.01% 

Proport. to signal level 

Linearity 

+0.01% 


Analog Subsystem 

Sampling 

Negligible 

Sampling rate & filter- 
ing adequate 

Multiplexer 

Transfer Accuracy 

+0.01% 

Proport. to signal level 

Settling Time 

= 0 


Rate 

*0 


Input Leakage 

= 0 


Sample and Hold 

Acquisition 

= 0 


Hold 

4 mV 


Aperature Delay 

= 0 


Offset 

= 0 

Corrected by software 

Gain 

*0 

Corrected by software 

Nonlinearity 

+0.01% 

Proport. to signal level 

Feedthrough 

+0.01% 

Proport. to signal level 
and frequency 

A/D Converter 

Resolution 

+1. 2mV 


Accuracy 

+1. 4mV 


Offset 

= 0 


Gain 

*0 



Table 4-9 SENSOR AND ANALOG SUBSYSTEM ERROR SUMMARY 




99 


4-6 SAMPLE ERROR ANALYSIS 


Orthogonality correction using the algorithmic method can 
be verified by computing expected error prior to correction and 
comparing measured system output with the error predicition. 
Assuming that the angle between the x and y sensors exceeds 90 
degrees- as in Fig. 4-13, we can proceed to compute error by 
noting the following relationships 

Hx = Hh Cos (+tjj) 
riy = Hh Sin (i!>) 

Hx i - Hx Cos e - Hy Sin e 

True Yaw = <,T = Cos' 1 
Computed yaw 

HxCose - HySine 

(HxCose - HySine) 2 + Hy 2 ) 2 

HhCosftCose - HhSinifrSine 

[(HhCosiJ/Cose - HhSinif»Sine) 2 + H 2 hSin 2 ipJ 2 

Using small angle approximations with e = 0.79° 

Cose = 1 and Sine * 0.014 


if>m = Cos 


-1 


= Cos 


-1 


then 


-1 


ij/m = Cos 


Computed error 


Hh CosiJ; - 0.014 Sinit 


[Hh 2 ( CosiJj — 0.014 Sini/;) 2 + Hh 2 Sin 2 <f/]^ 


Error = i^m - <p 


( 
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N 



Fxg . 4-13 ANGLE (X - Y)>90° 


We can now evaluate computed yaw angle (tj>m) given a part- 
icular yaw (ip) and the horizontal field vector {Hh) . Heading 


error. for horizontal field vector 

of 730 unirs at 

various yaw 

angles with pitch 

and roll angles 

of zero degrees 

is tabulated 

m Table 4-10 and 

plotted along with actual measured yaw error 

{data taken during 

1 experimentation 

. of Chapter V) 

m Fig. 4-14. 

• 

Computed 


Computed 

• Heading 

Error 

Heading 

Error 

(Degrees) 

(Degrees) 

(Degrees) 

(Degrees) 

90 

0.8 

290 

0.7 

70 

0.7 

270 

0.8 

50 

0.5 

250 

0.7 

30 

0.2 

230 

0.5 

10 

0.0 

210 

0.2 

350 

0.0 

190 

0.0 

330 

0.2 

170 

0.0 

310 

0.5 

150 

0.2 



130 

0.5 


Table 4-10 COMPUTED HEADING ERROR WITH Hh 


730 UNITS 




i 


COMPARISON OF COMPUTED AND MEASURED ORTHOGONALITY INDUCED ERROR 


Fo.g . 4-14 
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4-7 CONCLUSIONS 


The preceeding error analysis has identified potential 
error sources along with relative magnitudes of error to be 
expected. Magnetometer sensor and analog subsystem errors 
were identified and analyzed individually. During thus analy- 
sis it became apparent that errors due to sensor offset and 
nonorthogonality dominated and would severely limit total in- 
strument performance. The relative magnitudes of these errors 
and their mode of contribution would have degraded system cap- 
acity. 

By carefully characterizing the offset and orthogonality 
error it was determined that these systemmavic errors could be 
reduced by appropriate programming. A need to identify the 
extent of each error unique to the laboratory instrument im- 
posed a need to evaluate the instrument empirically. using 
earth's magnetic field and the laboratory test fixture (de- 
scribed in Chapter V) to provide control inputs each of the 
parameters was identified and measured. An algorithm with the 
empirically determined correction coefficients was included in 
the final system to reduce the error and to improve final system 
performance. The remaining potential error sources were tab- 
ulated and relative magnitudes noted. 

Processing errors due to register precision and truncation 
were analyzed by considering pertinent subroutines individually. 
It was noted that the relative error bounds add when multiply- 
ing variables with relative error. In addition, it was noted 
that error accrued during processing is proportional to sensor 
signal levels involved. The final uncertainty is then propor- 
tional to actual aircraft attitude with error increasing as 
displacement from level flight occurs. Computational error is 
also noted to increase at particular headings causing the error 
function to peak at specific yaw angles. 
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The sample error analysis clearly shows that a correlation 
between sensor nonorthogonality induced error and measured (un- 
corrected) data exists. By predicting and computing an error 
function prior to experimentally verifying the result we gain 
confidence that the sensor characteristics derived empirically 
m previous sections are correct. 
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CHAPTER V 


LABORATORY EVALUATION OF THE ATTITUDE INDEPENDENT 
REMOTE MAGNETIC INDICATOR AND HEADING INSTRUMENT 


5-1 INTRODUCTION 


• This chapter addresses laboratory evaluation of the micro- 
processor based computer designed to implement the heading 
measurement instrument. An integral part of this instrument 
was the three axis fluxgate magnetometer used to implement the 
attitude independent remote magnetic indicator of Chapter II. 

The laboratory evaluation was designed to investigate empiric- 
ally the effects of physical parameters that would otherwise 

» 

be impossible to assess. 

Although phenomena such as noise, magnetic field gradient, 
sensor orthogonality errors and offset errors can be predicated, 
combined effects on the proposed instrument and remote magnetic 
indicator are best evaluated in the laboratory. In addition, 
it was noted that errors due to sensor offset and nonorthogon- 
ality could be corrected by software included with the sample 
subroutine. Determination of the effectiveness of this cor- 
rection technique necessitated laboratory measurements of the 
errors (to determine correction constants) and comparison of 
data prior to and following corrections . 

The chapter begins by discussing laboratory test apparatus 
designed to evaluate the instrument. Actual data measured and 
recorded during experimentation is then presented in both tab- 
ular and graphic form to facilitate comparison and evaluation. 
Finally, the laboratory data is discussed and it is concluded 
that the remote magnetic indicator used with the heading mea- 
surement instrument results in a viable alternative to conven- 
tional heading measurement systems. The microprocessor based 
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computer implentation of the instrument has added unique sensor 
measurement correction ability that enhances performance of 
otherwise marginal sensors. In this manner limitations m 
systems performance that no w exist due to sensor inadequacy 
can be minimized without incurring the burden of using more 
expensive sensors. 

5-2 TEST APPARATUS 

A) Electronic Subsystem 

The microprocessor based computer (illustrated in photos 
5-1 and 5-2) was constructed on printed circuit boards consist- 
ing of a central processing card, two memory cards (2K bytes 
capacity each) and an output board. A separate analog subsystem 
card contained the multiplexer, jsample and hold, analog to 
digital converter and trimming potentiometers. The circuit 
cards were all organized with edge connectors and mounted ver- 
tically into a hand wired backplane assembly as shown in photos 
5-1 and 5-2. 

The card m the left foreground of photo 5-1 served as the 
output display with three seven-segment displays displaying 
significant figures of system heading. A small printed circuit 
in the right foreground of photo 5-1 contained potentiometers 
used to generate analog signals proportional to roll and pitch 
signals (simulating gyroscope outputs) . Cards shown vertically 
mounred in photo 5-2 can be identified from right to left as 
the analog subsystem, two memory cards and the central pro- 
cessing card. The large integrated circuit shown on the CPU 
card is the Signetics 2640 microprocessor. 

B) Sensor Assembly 


To evaluate the effects of combined aircraft pitch, roll 
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Photo 5-1 . MICROPROCESSOR BASED HEADING COMPUTER 



Photo 5-2 


CENTRAL PROCESSOR, MEMORY AND 
ANALOG SUBSYSTEM 
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and yaw a three axis gimbal apparatus was required. In add- 
ition, since angular measurements were required, a means of 
measuring angular rotation in each of the three exes was pro- 
vided. The gimbal apparatus as illustrated in photos 5-3 
and 5-4 was fitted with large protractors centered on the 
rotation axes. Pointers were provided to enable angular rot- 
ation measurements on the respective protractor scales. Since 
the angular precision on each protractor scale resolved angular 
position to 0.5 degrees, angular measurements to a resolution 
of at least 0.5 degrees were possible. Angular position was 
measured by estimating the decimal place of each measurement 
with accuracy to +0.5 degrees ensured. 


Since the three axis magnetometer (housed in the rectang- 
ular block of photos 5-3 and 5-4) measured ambient magnetic 
fields the test apparatus was constructed of nonferrous mat- 
This ensured that local fields due to residual mag- 
netic fields in the test apparatus would be minimized. In 
addition, since the material had low permeability, there would 
be little deformation of the local field causing error due to 
changing field gradient. 

The sensor package shown in photos 5-3 and 5-4 was physic- 
ally mounted such that the sensors were centered as close to 
the center of the gimbal as possible. This precaution ensur- 
ed that measurement error due to sensor translation was minim- 
lzed 1 . During instrument evaluation, the entire gimbal assembly 
and sensor were leveled and mounted in a Helmholtz coil assembly 
as illustrated in photo 5-5. Although the coils were not acti- 
vated during experimentation, the rotations in heading were 

Since the local magnetic field has a nonzero gradient, field 
measurements include a component due to translation of the 
sensor axes. This component of measurement produces unaccept- 
able error in a system designed to measure field components 
that change due to rotation. 
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Photo 5-3 


MAGNETOMETER SENSOR MOUNTED ON 
GIMBALLED TEST FIXTURE 



Photo 5-4 


SENSOR AND GIMBAL ASSEMBLY WITH PROTRACTORS 
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carefully controlled since the gimbal assembly was an integral 
part of the Helmholtz coil fixture with the vertical rotation 
axis serving as the system yaw axis. 
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Photo 5-5 TEST FIXTURE MOUNTED IN HELMHOLTZ COIL ASSEMBLY 


5-3 HEADING MEASUREMENTS WITH NO OFFSET CORRECTION 

By maintaining heading of the text fixture constant (no 
rotation about the vertical axis) and varying both pitch and 
roll angle, the instrument display was observed to vary. This 
variation gave a direct measure of instrument error since a 
constant heading was maintained and a constant display was to 
be expected. 

Data variations were recorded in Tables 5-1 and 5-2 and 
plotted on Figures 5-1 and 5-2. With only +10 degree variation 
in pitch combined with +30 degree variation in roll we note 
that the heading display varies 14 degrees. Obviously, 


Heading Computed (Degrees) 



Pig. 5-1 


Heading Computed (Degrees) 



Roll Angle (Degrees) 


Fig. 5-2 HEADING COMPUTED AT A FIXED YAW ANGLE WITH VARYING 
G AND 4 (NO OFFSET CORRECTION) 
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instrument operation indicated excessive error requiring more 
elaborate sensors or correction of a sensor inadequacy . 

5-4 HEADING MEASUREMENTS TO INVESTIGATE ORTHOGONALITY ERROR 

System performance was evaluated by initially aligning the 
sensors with zero pinch and roll angle. Sensor Z was position- 
ed vertically with positive direction downwards. By observing 
the Z axis output 2 as the test fixture was rotated about the 
vertical axis, adjustments were made m pitch and roll angle 
to minimize coning of the Z axis. Angular measurements on 
the respective roll and pitch axis protractors were then made 
to establish the initial reference attitude angles. 

Heading measurement accuracy was evaluated by rotating 
the test fixture in the horizontal plane until the display 
flickered between (XX9) and (XX9+1) . The rotation was then 
continued a very small amount until a steady display (multiple 
of 10 degrees) was -observed 3 . Measurements ranging from 0 to 
350 degrees were made by recording angular position required 
to produce specific heading data displays. Sets of data were 
recorded at various combinations of pitch and roll then tab- 
ulated in Tables 5-3 through 5-8. Relative error was computed 
by determining angular position expected at each display value 
and then computing the difference m angular positions. Errors 
at the roll extremes of + 44 degrees are plotted for pitch 
angles of plus and minus 20 degrees on Fig. 5-3 through 5-6 
inclusive. 


2 A special subroutine was. used to display Z axis data directly 
in BCD format on the seven bar output display. 

3 This measurement technique ensured that all heading measure- 
ments were made identically. In addition, error due to. system 
imprecision was reduced. 
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Data in Tables 5-5, 5-6 and. Fig. 5-3, 5-4 were recorded 
with no sensor orthogonality error correction implemented. Data 
in Tables 5-7, 5-8 and Fig. 5-5, 5-6 was recorded with the 
sensor orthogonality correction implemented. Comparison of 
these data indicate that considerable improvement in accuracy 
is achieved by correcting sensor orthogonality error. 

5-5 CONCLUSIONS 

Laboratory evaluation of the heading measurement instru- 
ment has shown that the algorithms developed in previous 
chapters are viable. Operation of the device in a laboratory 
environment has enabled empirical evaluation of the system 
under adverse combinations of noise, field gradient and sensor 
plus instrument error sources. 

Test apparatus described in section 5-2 served to enable 
controlled simulation c5f roll, pitch and yaw rotations. The 
apparatus was nonmagnetic m nature and contributed insignifi- 
cant error due to field pertebation . Mounting of protractors 
and pointers on the test apparatus made angular measurements 
possible to a precision of at least +0.5 degrees. 

Effects of sensor offsets were evaluated in section 5-3 by 
recording system heading computations when only roll and pitch 
varied. 'Since the variations in Figures 5-1 and 5-2 prior to 
offset correction exceed the maximum excursions of Figures 5-3 
and 5-4 by at least a factor of two (angular excursions m 
first set also less than in the record) and we note that offset 
errors were corrected prior to recording data m the second set 
of‘ data, we conclude that offset m magnetometers can be a 
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Heading 

Displayed 

(Degrees) 


10 

30 

50 

70 

90 

130 

150 

170 

190 

210 

230 

250 

270 

290 

310 

330 

350 


PITCH ANGLE 0 DEGREES 
ROLL ANGLE 0 DEGREES 

Angular 

Position 

(Degrees) 

275.3 

295.2 

315.4 

335. 5 

355.5 
34.8 
54.3 

74.0 

94.0 

114.4 
135.0 

155.5 

174.7 

195.5 

215.2 

235.2 

254 . 7 


Relative 

Error 

(Degrees) 

0.3 

0.2 

0.4 

0.5 

0.5 

- 0-2 

-0.7 

- 1.0 


- 1 . 0 . 


- 0.6 

0.0 

0.5 

-0.3 

0.5 

0.2 

0.2 

-0.3 


Table 5-3 REFERENCE DATA MEASUREMENTS OF HEADING 
TAKEN WITH NO ORTHOGONALITY CORRECTION 
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PITCH ANGLE 0 DEGREES 
ROLL ANGLE 0 DEGREES 


Heading 

Angular 

Relative 

Displayed 

Position 

Error 

, (Degrees) 

(Degrees) 

(Degrees) 

10 

276. 6 

-0.4 

30 

296.6 

-0.4 

40 

306.7 

-0.3 

50 

316.7 

-0.3. 

60 

327.2 

+ 0.2 

70 

337.0 

0.0 

90 

355.9 

-1.1 

130 

37.0 

0.0 

160 

67.2 

0.2 

190 

96.8 

-0.2 

220 

126.9 

-0.1 

250 

157.0 

0.0 

230 

186.9 

-0.1 

310 

216.6 

-0.4 

340 

246.9 

-0.1 

350 

256.4 

-0.6 


Table 5-4 HEADING MEASUREMENTS WITH OFFSET AND 
ORTHOGONALITY CORRECTIONS MADE 



Roll = 44° 


Roll = 20° 


Roll = -20° 


Roll = -44° 


f Heading 
* Displayed Angular Angular 

(Degrees) Position Error Position 


10 

276.3 

1.3 

276.5 

30 

296.3 

1.3 


40 



306.1 

50 




60 

325.8 

0.8 


70 



335.6 

90 

355.0 

0.0 

354.0 

130 

34.3 

0.7 

34.3 

160 

64.5 

-0.5 

64.0 

190 

94.3 

-0.7 

94.6 

220 

125.2 

0.2 

125.3 

250 

156.6 

1.6 

156.2 

280 

186.6 

1.6 

186.8 

310 

216.8 

1.8 

216.8 

340 

246.0 

1.0 

246.8 

350 





Error 

Angular 

Position 

Error 

Angular 

Position 

Error 

1.5 

275.0 

0.0 

274.0 

-1.0 




293.6 

-1.4 

1.1 

304.0 

-1.0 

313.8 

-1.2 

1.6 

334.3 

-0.7 

333.8 

-1.2 

-1.0 

353.1 

-1.9 

353.0 

-2.0 

-0.7 

34.6 

-0.4 

34.4 

-0.6 

1 

W 

♦ 

O 

65.0 

0.0 

65.0 

0.0 

1 

o 

• 

95.0 

0.0 

95.5 

0.5 

0.3 

126.0 

1.0 

126.5 

1.5 

1.2 

156.7 

1.7 

157.0 

2.0 

1.8 

186.5 

1.5 



CO 

• 

r-4 

215.9 

0.9 

215.9 

0.9 

1.8 

245.0 

0.0 

224.9 

-0.1 


TABLE 5-5 HEADING MEASUREMENTS AT PITCH - 20° WITH NO 
ORTHOGONALITY CORRECTION 
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Fig „ 5-3 HEADING MEASUREMENT ERROR AT PITCH = 20° 
(NO ORTHOGONALITY CORRECTION) 
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Heading 

Roll = 

44° 

Roll = 

20° 

• Roll = 

-20° 

Roll = 

-44° 

Displayed 

Angular 


Angular 


Angular 


Angular 


(Degrees) 

Position 

Error 

Position 

Error 

Position 

Error 

Position 

Error 

10 

276.1 

1.1 

275.1 

0.1 

274.2 

-0.8 

273.7 

-1.3 

30 

296.2 

1.2 

295.9 

0.9 



293.6 

-1.4 

50 

316.4 

1.4 

316.1 

1.1 



314.1 

-0.9 

70 

336.5 

1.5 

336.4 

1.4 

335.5 

0.5 

334.6 

-0.4 

90 

355.5 

- 0.5 

355.3 

0.3 

355.0 

0.0 

354.2 

-0.8 

130 

36.0 

1.0 

36.4 

1.4 

135.7 

0.7 

135.7 

0.7 

150 

55.3 

0.3 

55.9 

0.9 



55.8 

0.8 

170 

74.5 

-0.5 

75.6 

0.6 



75.5 

0.5 

190 

94.9 

-0.1 

95.3 

0.3 

95.5 

0.5 

95.5 

0.5 

210 

115.0 

0.0 

115.4 

0.4 



115.6 

0.6 

220 





125.5 

0.5 



230 

135.0 

0.0 

135. 6 

0.6 



136.0 

1.0 

250 

155.6 

0.6 

155.4 

0.4 

155.5 

0.5 

155.7 

0.7 

270 

174.8 

-0.2 

174 . 8 

-0.2 



175.0 

0.0 

290 

195.6 

0.6 

195.4 

0.4 



195.3 

0.3 

310 

215.5 

0.5 

215.3 

0.3 

214.6 

-0.4 

214.5 

-0.5 » 

330 

235.7 

0.7 

235.1 

0.1 



234.2 

-0.8 * 

350 

255.6 

0.6 

255.2 

0.2 



253.8 

-Is 2 


Table 5-6 HEADING MEASUREMENTS AT PITCH = -20° 
WITH NO ORTHOGONALITY CORRECTION 
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ffl 



fxFm 


[•left] 



UREMENT ERROR AT PITCH = -20° 
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Heading 

Displayed 

(Degrees) 

Roll = 

Angular 

Position 

0° 

Error 

Roll = 

Angular 

Position 

44° 

Error 

Roll = 

Angular 

Position 

-44° 

Error 

20 

287.3 

0.3 

287.5 

0.3 

286.3 

-0.7 

40 

307.2 

0.2 

307.5 

0.3 

306.0 

-1.0 

60 

327.0 

0.0 

327.9 

0.7 

325.9 

-1.1 

80 

346.3 

-0.7 

347.5 

0.3 

345.6 

-1,4 

90 

355.9 

-1.1 

357.0 

0.2 

356.2 

-1.8 

140 

46.9 

-0.1 

47.0 

-0.2 

46.5 

-0.5 

160 

67 . 0 

0.0 

67.3 

0.1 

67.1 

0.1 

180 

86.6 

-0.4 

86.5 

-0 . 7 

87.0 

0.0 

200 

107.0 

0.0 

107.3 

0.1 

108.1 

1.1 

220 

127.6 

0.6 

127.3 

0.1 

128.3 

• 1.2 

240 

147.2 

0.2 

147.4 

0.2 

148.5 

1.5 

260 

167.0 

0.0 

167.1 

-0.1 

168.5 

1.5 

280 

187.3 

0.3 

187 . 1 

-0.1 

188.5 

1.5 

300 

207.2 

0.2 

207.1 

-0.1 

208.0 

1.0 

320 

227.3 

0.3 

227.0 

-0.2 

227.5 

0.5 

340 

247.4 

0.4 

247.3 

+0.1 

247.4 

0.4 

0 

266.4 

-0.6 

267.0 

-0.2 

266.0 

-1.0 


Table 5-7 HEADING MEASUREMENTS AT PITCH =20° 

WITH OFFSET AND ORTHOGONALITY CORRECTION 
MADE 




-2. 0 _L 

Fig. 5-5 HEADING MEASUREMENT ERROR AT PITCH - 20° 

(OFFSET AND ORTHOGONALITY ERROR CORRECTED) 
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122 


123 


Heading 

Displayed 

(Degrees) 


Roll = 44° 
Angular 

Position Error 


Roll = -44° 
Angular 

Position Error 


20 

287.0 

-0.4 

286.5 

-0.3 

40 

307.0 

-0.4 

306.7 

-0.1 

60 

327.0 

-0.4 

326.8 

0.0 

80 

347.0 

-0.4 

346.8 

0.0 

90 

356.7 

-0.7 

355.6 

-0.2 

140 

47.8 

0.4 

48.0 

1.2 

160 

67.7 

0.3 

68.0 

1.2 

180 

87.0 

-0.4 

87.8 

1.0 

200 

107.6 

0.2 

108.0 

1.2 

220 

127.7 

0.3 

12 8.0 

1.2 

240 

147.8 

0.4 

147.4 

0.6 

260 

167.6 

0.2 

167.2 

0.4 

280 

187 . 6 

0.2 

187.0 

0.2 

300 

207.1 

-0.3 

206.5 

-0.3 

320 

227.1 

-0.3 

226.6 

-0.2 

340 

247.1 

-0.3 

246.9 

0.1 

0 

266,3 

-0.9 

266.0 

-0.8 


TABLE 5-8 HEADING MEASUREMENTS AT PITCH = -20° WITH 
OFFSET AND ORTHOGONALITY ERROR CORRECTED 
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major error source 4 . Additionally, we note that the correction 
'of offset error in sensors has been successful. Experimental 
results have verified that not only can offset errors be deter- 
minded (Chapter IV) , but a suitable algorithm can be implement- 
■ed in the computer to improve system operation. It is postu- 
lated that offset error correction can be extended to include 
correction of varying offset values (functions of temperature 
and supply voltage) by monitoring error causing' variables 
(example temperature) and computing correction constants prior 
to offset correction as above. 

Errors induced by sensor nonorthogonality were predicted 
in Chapter IVsection 4-2 and verified by plotting expected 
error along with measured error in Fig. 4-14. The curves of 
Fig. 4-14 were plotted for heading rotations with no pitch or 
roll angle. To evaluate system performance and the effect of 
orthogonality error with combined- angular rotations , measure- 
ments of heading error were plotted m Fig. 5-3 thorugh 5-6 
inclusive. 

Comparison of these data indicate that maximum excursions 
of error as a function of heading are significantly less when 
orthogonality corrections are made. It is also postulated that 
data could be improved further by similarly correcting ortho- 
gonality error m te Z axis sensor 5 . 

In summary, the experimental evaluation has provided in- 
sight into the operation of an attitude independent remove mag- 
netic indicator and heading computer in the "real world" 


4 This corroborates the observations predicted during error 
analysis in Chapter IV. 

5 We note that the error excursions are functions of pitch 
and roll and that Z axis data is used in the rotation algorithm. 



126 


environment complete with all contributing error sources. The 

error analysis evolved during development of the system has 

proven adequate in that an operational system was developed. 

Major error sources were measurable as predicted and the means 

* 

of reducing their effects were successfully implemented. Cor- 
rection of sensor offset and orthogonality error required an 
empirical evaluation of the respective sensor. These evalua- 
tions were performed, the errors characterized, correction co- 
efficients determined, and correction algor itms implemented.- 

Successful implementation of these corrections was evi- 
denced by significant reductions m system error. The correct- 
ion methods presented can be extended in future with the net 
result that less demand is required of physical sensors if the 
sensor parameters can be established empirically prior to 
completion of instrument design. Utilization of a microproces- 
sor in the instrument has added the computational flexibility 
required to facilitate accommodation of sensors with varying 
error magnitudes . 



APPENDIX A 


This appendix lists the instruction set of the Signetics 
2650 microprocessor chip used to implement the heading instru- 
ment. 
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loao/store instructions 

Unylh (bytes) 

BJRVr (+)a 

Brand) on Incrementing Register Absolute 

3 

LOD Z 

r 

Load Register Zero 

1 

BDltR,r (*)u 

Bunch on Decrementing Register Relative 

o 

LOIR r 

V 

Load Immediate 

2 

BDRA,r {+)a 

Brandt on Decrementing Register M>solute 

3 

LODlt.r 

l*)a 

Loud Relative 

2 

BXA (*MU) 

Bunch Indexed Absolute Cncond mortal 

3 

LOOA,r 

(*>a(.X) 

I oad Absolute 

3 

ZBRR <*ja 

Zero Brandi Relative, Unconditional 

2 

STR2 

r 

SLOre llegibLer 2u.o 

1 




STRR.r 

(*)u 

Store Relative 

2 

SUBROUTINE BRAN CH/flE TURN INSTRUCTIONS 


str A,r 


Store ALfooluio 

3 

BSTR.v (*p 

Branch to Subroutine on Condition 

2 






True, Relative 


arithmetic instructions 


BSfR.v (*)a 

Branch to Subroutine on Condition 

2 

AD 02 

r 

Add to Register 2eio 

1 


FuLe, Relative 


ADOJ r 

V 

Add Immediate 

2 

BSTA.v (+)u 

Branch to Subroutine on Condition 

3 

ADDR,r 

(* M 

Add Relative 

2 


True absolute 


ADD V,r 

(*ja(,X) 

Add Ab>[)tuii 

3 

USFA t v (*)a 

Branch to Subroutine on Condition 

3 

SUB/ 

r 

bubtrai i from Rcgibter Zero 

1 


False Absolute 


subi r 

V 

Subtract immediate 

2 

BSN R,r ( * )u 

Branch to Subroutine on Non Zero 

2 

SUBR.r 

(*)a 

Subtract Uelaiive 

2 


Register, Relative 


bUBA.r 

( + )a(.X) 

Subtract Absolute 

3 

BSNA,r (*)u 

Brunch to Subroutine on Non Zero 

3 






Register, Absolute 


LOGICAL INSTRUCTIONS 


BSXA (*)*(,*) 

Brandi to Subroutine Indexed Unconditional 

3 

ANDZ 

r 

And lo Register Zuo 

1 

IlErC.v 

Return From Subroutine Conditional 

1 

AND! r 

V 

And Immediate 

2 

UCTC,v 

Return From Subi online and Fntible 

1 

ANDlljr 

<*)•* 

And Relative 

2 


Interrupt, Conditional 


AND.A.r 

(*)aLX) 

And Absolute 

a 

ZBSR <*),U 

2t.ro Branch to Subiouiinfr 

2 

10 IU 

r 

Inclusive or to Register Zero 

i 


Relative, Unconditional 


10 III, r 

V 

Inclusive or Immediate 

2 




10RR,r 


Indusive ot Relative 

2 

PROGRAM STATUS INSTRUCTIONS 


IOU \,r 

(+)a( > X J 

Inclusive or Absolute 

3 

LPfaU 

Load Progi*u« Status, Upper 

1 

LOR L 

r 

Exclusive or to Register Zero 

l 

LPSL 

load Program Status, Lov^r 

1 

liOKI r 

V 

Exclusive or Immediate 

2 

SPSU 

Score Program Status, Upper 

1 

LOKJt.r 


Exclusive or Relative 

2 

SPbL 

Stoie Program Status. Lower 

l 

UUKA,r 


CxcJuMve or Absolute 

3 

CPSU V 

Cltai Program Status, l ppei Selective 

2 





CPSl v 

CLur Program Status, Loner Selective 

2 

COMPARISON INSTRUCTIONS 


PPSU V 

Preset Program Status, Upper Solectne 

2 

COM2 

r 

Compare to Register Zero 

1 

PPSL v 

Prcbd Program Statu*. Lower Selective 

2 

COMI r 

V 

Compare Immediate 

2 

TPSU V 

Test Program Status Upper selective 

* o 

COMICr 

*^|U 

i oiupdre Urdative 

2 

TIM- v 

'kst Pi o gram Status J owtf Selective 

2 

COM Mr 

Xj 

Compare Absolute 

3 








INPUT/OJTPUT INSTRUCTIONS 


ROTATE INSTRUCTIONS 

Length Ibytut) 

WKTD.r 

Write Data 

l 

ItUK.r 


Route Rcgibter Right 

l 

RKDD.r 

JU.id Data 

l 

It It L,r 


Route Rcgibter Lett 

1 

WltTC.r 

Write Control 

l 





RCDG'.r 

Read Control 

1 

BRANCH INSTRUCTIONS 


WftTE.r v 

Write LxtvmUd 

2 

B( lR,v 


Brand \ on Condition 'lYuo Relative 

2 

RUDix,r v 

Read Extended 

0 

aci n v 

(*;a 

Bmmh on Condiiion False Relative 

2 




BC 1 A,v 

(♦)a 

Brain h on Condition 'IVue Absolute 

d 

MISCELLANEOUS INSTRUCTIONS 


}IC f \ v 

(*hl 

Brand) on Condition False Absolute 

3 

HAJ T 

Hall Enter Watt Stab 

1 

BKNK,r 

(*|a 

Brandi on Register Non-Zero Relative 

l 

DA 11 r 

IX'cmuI AdjUbt Register 

1 

URN V,r 

( 

Br.uu h on Register Non Zero Absolute 

3 

TM1 r v 

'D’bL Under M*uk 1m mediate 

2 

BUUCr 

(*)a 

Brandi on Incrementing Register Relative 

o 

NO? 

No Operation 

l 


H* 

tv) 

OO 



REPRODUCIBILITY OP THE] 
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APPENDIX B 


This appendix contains alisting of the assembly language 
program used to implement the remote magnetic indicator head- 
ing algorithm. The program was assembled on ^ the A2650 cross 
assembler program operational on the HP 2100 computer at the 
University of Santa Clara. 
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LINE 

J 

4 

5 

6 
; 
a 

9 

10 

n 

ia 

13 

14 

15 

16 
17 
ia 

19 

20 
ai 
22 

23 

24 

25 

26 
dl 

aa 

29 

jo 

ji 

32 

33 

34 

35 

36 

37 
Ju 
jy 
no 
M 
*♦2 
43 

45 
*6 
*./ 
4 a 
4V 

50 

51 
a? 
bj 

b4 


HIP ASSfMn|_. H VT WS 1 f ‘N SfU LI VLI 1 HI- AO [ NG 1 N<> I RUMl- NT ASSEMLY PROGRAM l‘J7h PARf 2 


AOOP L XhL HI LV h 3 <U 4 t riuGR SOURCE 


0000 

MO 

EOU 

0 

rifnu 

W1 

LOU 

1 

0 u 02 

HZ 

EUU 

2 

OunJ 

RJ 

Eou 

J 

01)0 0 

ON 

t OU 

0 

on«o 

L 

EDO 

0 

0001 

P 

LOU 

1 

<u>n2 

N 

EDO 

2 

do no 

f-Q 

LOU 

0 

0 0 1) 1 

GT 

EuO 

1 

(1001 

r 

LOU 

l 

00()2 

LT 

EOU 

2 

0003 

UN 

Eou 

3 

04)02 

COM 

EuU 

H*02» 

no i 0 

MS 

Lou 

H 1 1 0 ♦ 

o^na 

wC 

Li HI 

N Hjrt * 

00 or 

LUWR 

LOU 

H* OF » 

oofo 

UPPR 

EOU 

H * F 0 1 

0040 

FLAG 

Eou 

H ) 4 0 ' 

00P0 

II 

Enu 

Hi20« 

OUFf 

NSAM 

EuU 

H*FF * 

0001 

MSB 

tOu 

1 

0002 

LSB 

LOU 

2 

0004 

THTA 

LOU 

4 

ODnB 

PHI 

too 

8 

onto 

Sri I T 

LOU 

16 






data 

STORAGE IN HAH 



ORG 

H»450* 

04SO 

DATA 

RES 

1 n 

44SA 

OATH 

RES 

in 

04*4 

LUPS 

RES 

1 

04nb 

TEMP 

HFS 

2 

i)4*>7 

SIND 

HLb 

2 

04*9 

SINR 

RES 

2 

i)4*>b 

hn2 

RtS 

2 

046U 

SIGN 

RtS 

j 

04s£ 

nX? 

RES 

2 

04 /0 

HY2 

RES 

2 

04 72 

nid 

RES 

2 

44 74 

UUAO 

RES 

l 

4475 

COSR 

RES 

2 

44 77 

RSLT 

RtS 

2 

44/9 

OPR] 

RtS . 

2 

U / o 

OPR2 

RLS 

2 

44 /U 

STEM 

RtS 

l 

04 /L 

SUES 

RFS 

1 

04 7f 

DA T X 

RFS 

2 

04H1 

Oaty 

RtS 

2 

04 13 

I’HU 

RtS 

2 

04 Ab 

COST 

Rf S 

1 




OUTPUf 

OUTPUT 


«<*#» 

RAW Data 
CORRECTED DATA 


HM SUUAh£u J AT A 
riA bOilAPfiJ DATA 


OJ 
» I — 1 


HOO<r SI 30 VJ 3VNI0IH0 

30 ^imsoiigog^ag 
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LiNF, 

AOUP 

1 A *L 

01 

H2 

*3 

tl4 F hPOH 

<■ ounce 



bS 


/Wh6 




- 

SIGx 

RcS 

1 

bfe 


0‘«h7 





S1GY 

kts 

1 

b/ 


04AP 





SIG^ 

MfcS 

1 

b« 


0**.9 





AHUF 

HF’b 

l 

b9 


f»**H A 





btJUF 

wts 

1 

60 







CtU)F 

PFb 

l 

6) 


n^MC 





DHllF 

RF S 

1 

02 


H4-*0 





THUF 

HFS 


6 ) 







ObET 

RES 

6 

C>4 


*W >4 





sosr 

HI- S 

6 

6b 








66 







uq»H 

CuNbT an ( s 

wQ M ao»t» 

67 








ORG 

H • 5D0 f 

6a 

Oaurt 

ObOO 

00 

01 

03 

04 

SlN 

DATA a 

Hi 00*0 1 ,0 i*04t06*07»Q9,OA,OCiO£»OF ■ 




v6 

07 

69 

UA 







OL 

0^ 

OF 





69 

0 6 OH 


1! 

W 

1C 

Is 


PaTa 

Hil 1 *12. U, Ib* 1 7* 19.1A. 




U 

19 

1 A 

lr 







ID 

K 

til 





70 

0616 


it 

2<* 

a5 

ti 


pa r a 

hi2?.24.2S*2/.2H,2A*2ti,2f>,2E* JO * 3 l « 




/a 


2b 

2n 







1 1 

10 

Jl 





/ l 

ob<M 


*13 

J-. 

16 

Ih 


DAT** 

n<3 <i lb. 3s ♦ 36. ly i Jo* J( *11- , JF*A| ( 42» 




i9 

I 1 

1C 

IF 







IF 

Cl 

42 





i? 

Ob?C 


t>4 

4 s 

4/ 

sn 


PaFm 

H*44*4b»4/t4Hi4A*4H»4D*4F . bO » b l *b J * 




4A 

4 H 

4 1) 

4F 







so 

6| 

S3 





/ J 

0bl7 


l iA 

56 

6 7 

bo 


Data 

h i b4 «S6 *b / »59« bA *bC*50 1 61 Of 6* 




AA 

bF 

SI) 





/4 

Gb^F 


bf 

64 

61 

61 


D*t« 

h > 5F toO.bl *63.64*66. 6 / *o9,6a.6Ci6u* 




44 

6r» 

6/ 

69 







4 A 

S^ 

->U 




h.6F * /0 * 7 1 *73.74. fbJjf* 7H.7a,7H*7c» 

/b 

Ubs9 


6f 

?■> 

71 

/l 


DATA 




?4 

/S 

// 

7rt 







;a 

/* 

'C 





/ 6 

ObV* 


7E 

/f 

HO 

Hp 


DATA 

h • 7r » 7F * «o * A2 * H J *H4 * 36* 67 * faA * A A i 06 « oC • 




Hi 

a 4 

4 6 

H 7 







a a 

*u 

4M 

ar 




// 

OboO 


Hfc 

P 

S'J 

9/ 


DA ro 

H.bl- .sF .9,1*92.93. 94.96*97,9 i.yO.Vn' 




9j 

94 

*S 

9/ 

, 






OH 

69 

#b 





/a 

DbnH 


‘>C 

* 1 

4t 

9f 


data 

H.9C.9D.9F »9F . Al * A2.A 1* A4, A6, A7 .Ac » 




A 1 

\s 

A ) 

A4 







ah 

A 7 

0 ft 





79 

05/6 


Av 

A A 

AH 

Ai) 


DaTa 

H 1 1 1 * » 41 A . A* i * AO * At • A F #H'1**M ,d2iHl*pb* 




At 

P J 

*0 

*1 







»V 

-ll 

ib 





60 

obai 


»« 6 

rt7 

HH 

riM 


Data 

Hips, ,7,H4*p9*iiA»F)D*ac*«P*HF*C0.CI *C2* 




HA 

W 

MC 

U|) 


' 





,-(■ 

( A 

Cl 

C/ 




61 

obai> 


f i 

Co 

cs 

C6 


OA I A 

Hit Uf4.Cb.Cb»r7.CH*C9.t A ,C*1»CC * CiH 




C/ 

( H 

( M 

f A 







( H 

cc 

C i> 







U> 

to 
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UnE 

ADOR 

LAUL 

8 1 

B2 

83 

04 ERROR 

SOURCE 

82 

0S98 


ce 

CF 

00 

01 





D2 

03 

03 

D4 





DS 

06 

L) 7 



d3 

0bA3 


DB 

09 

1)9 

Da 





08 

DC 

OD 

DO 





OE 

DF 

EO 



04 

OSaE 


FI 

El 

C2 

El 





£3 

E4 

ES 

£6 





E6 

E7 

EB 



U5 

05B9 

flbM9 

FF 

FF 

FF 

FF 

COS 




FF 

FF 

FF 

ff 





FF 

FF 




06 

05C3 


FF 

FF 

FF 

ff 





FF 

FF 

T E 

ff 





FE 

FE 




07 

OSCD 


FE 

FD 

FO 

Fn 





FD 

FC 

FC 

Fc 





FC 

Fn 




UB 

05D7 


F0 

FB 

FB 

Fa 





FA 

FA 

F9 

F 9 





F9 

FA 




09 

ObEl 


ffa 

F7 

F7 

F7 





F 6 

F 6 

F5 

FS 





FA 

Fa 




90 

05EB 


F 4 

F3 

F3 

F? 





F2 

FI 

FI 

Fn 





EF 

EF 




91 

OSFS 


EE 

EF 

eo 

£d 





EC 

Eb 

ED 

Ea 





E A 

E9 

EB 

EP 


92 

0601 


£7 

E6 

E6 

E5 





£A 

E3 

E3 

E? 





El 

El 




93 

060M 


EO 

DF 

BE 

00 





DO 

OC 

UB 

Da 





09 

09 




94 

0615 


OB 

07 

06 

05 





04 

03 

03 

D? 





m 

On 




95 

Ob l F 


LF 

CF 

(0 

c r 





CD 

CA 

C9 

CA 





C7 

CA 




96 

0629 


C5 

C4 

C3 

c? 





Cl 

Co 

HP 

be 





DO 

DC 




97 

0633 


wd 

BA 

H9 

Bp 





0 7 

U6 

bS 

83 





H2 

ai 




90 

06 JD 


RO 

AF 

aE 

ad 





aB 

AA 

A9 

aa 





A7 

A* 

A A 



99 

06A fl 


A3 

A2 

A 1 

9f 





9E 

90 

9C 

90 



DATA H*C£»CF»DO*Dl » D2 *03 *03*04 » 05 *06*07 • 

DATA HiDb»C9*09*0A*QB*0C »DD iDO , DE *DF *£0 1 

OaT A H»E1 ♦El*£2tE3tE3*£4*E5*£6 l E6*E7*ea» 

DATA HifF.f F .FF,FF*FF*FF,FF*Ff ,FF*FF» 

DATA hiFF*Ff *FF *FF*FF *FE *FE*FE *F£ *FE ♦ 

DATA HiF£*FD*fD,FD.FO,FC.FC.FC*FC*FB» 

DATA HiFB»FB*FB*FAiFA*FA*F9*F9*F9»FB» 

DATA HiF0*F7*F7*F7iF6*F6*F5iF5*F4*F4» 

DATA h.F4*F3*F3*F2*F2*Fl*FWF0.£F*£F» 

DATA HtEE*FE*ED*E0*EC*E8*£B*£A*EA*E9*EB*£B 

DATA H*E7*E6*E6*E5*£4*E3*E3*£2*El*El • 

DATA rl t EO • DF tOF * 00 * OD *0C » Dri *DA » D9 ,09 * 

DATA Hi 08 *07*06, Ob *04 *03. 03 *02,01 *00* 

DATA H * Cf » C£ t CD •CC,C0*CA»C9*CA*C?,C6' 

DATA HiCb*C4*C3»C2*Cl p CO • OF . BE , 0D* OC « 

OATA m0ri*OA*H9»08*tl7*06*05iB3*8?*Hl • 

DATA H*B0*aF*AE,AD*AD,AA*A9*AB,A/,A6.A4« 

DATA ti t A J « A<£ * A 1 *9f*9£, 90*90*98*99,90,97, 95 
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LINE 

ADOR 

t Aril. 

Hi 

a? 

S3 

B 4 f RkUK 

SOURCE 



100 

0664 


94 

93 

9? 

90 


DATA 

Hi94.93.9? t 90»6F»8Ei8C»8FUBA,aBi 




Hf 

BE 

BC 

Bfl 







HA 

H,i 






101 

065E 


rt 7 

HA 

84 

83 


DATA 

H<87,H6.84*83*B2*80i7F ► ?F , 70*70* 




82 

BO 

7F 

7f 







7C 

7B 






102 

0668 


7 A 

70 

77 

75 


DATA 

Hi7a» 7H*/7*7S*74*73*71 « 7o,6E,6D*6C* 




74 

71 

71 

7 0 







6t 

60 

6C 





103 

06 71 


6A 

69 

67 

66 


OAT a 

H * 6 A * t>9 ,6/i66«64,63t6l * 60 * bF , SO * 5C * 




64 

63 

61 

60 







5F 

SO 

bC 





104 

06 7t 


b A 

59 

57 

Sa 


DATA 

H,54,59,57»56,S4*53*51*5O,4F,40» 




54 

S3 

51 

50 







4b 

40 






10b 

06H8 


4rl 

4 A 

4 a 

47 


DAM 

h , 4H * 4A » 4H 1 47 *45*44,4? *4 1 » Jf , 3E , 3C * 




45 

44 

42 

4 l 







IF 

3F 

3C 





106 

0691 


3B 

J9 

38 

3 A 


DAT w 

H » 38 . iy * 3* • 16 . 35 . 3 3 ♦ 3 1 . 3 0 * 2£ , 20 * 




35 

33 

Ji 

30 







?t 

20 






107 

U69f> 


2b 

2A 

?a 

27 


0 A 1 A 

H *28. PA, ?Hi 27 ,25*24,22*20, IF, ID* 




2b 

24 

22 

20 







IF 

10 






108 

0647 


1C 

1A 

19 

17 


DATA 

rt»lC*tA.l9, 1 7. IS# 14,12*11 ,0F»OE* 




15 

l 4 

12 

11 







OF 

Of 






iu9 

(J6r)l 


oc 

Oa 

09 

07 


DATA 

N* OC* OA*09 * 07 » 06 * 04* 0 1 * 01 « Ofl » 




06 

0/* 

03 

01 







00 







no 

06HA 

D6RA 

7F 

7H 

77 

72 

COSM 

DATA 

H,7F *7H* 77*72 *6E *69,65.61 ,5c*Sfl»b4* 




6b 

69 

65 

6 1 







bC 

5H 

64 





m 

06CS 


50 

4‘, 

4 7 

43 


DATA 

H*50 ,40,47* 43* 3F « JC ► 3H# J4 . 31 *2D» 2A * 




3F 

JC 

38 

34 







31 

20 

2A 





112 

o6no 


27 

21 

20 

It) 


data 

H*2/*?3.2o, ID. lfl*lQ*lS*13,U ,O£,0L» 




lb 

1« 

15 

11 



* 




11 

OF 

OC 





113 

06i)R 


OB 

09 

07 

Oa 


DATA 

H» OH* 09,0 7 i 06 » 04 * 03* 02*01, 01 *00*00 1 00* 




04 

03 

02 

01 







01 

0 0 

00 

00 







00 







1 1 4 

06E8 

OOF B 

FF 

OH 

1? 

9F 

COSL 

DATA 

H*FF* 88.12,9b »2E*C5* 63 *08,87* 72.38.0C* 




2 F 

C5 

63 

Oa 







R7 

7? 

38 

or 




lib 

U6F 4 


\ 0 

El 

F 8 

*F 


OA f A 

mF0*tJ.LB.FF*2ti*t,C,C3*J2,B9,59, ib,tC* 




2B 

6C 

C3 

3? 







H9 

59 

15 

£C 




116 

0700 


FO 

F? 

22 

f? 


DATA 

M«EO,f 2 *22*72 *E2* 73,26 *FB,F4. 10,52 *dB» 




f2 

73 

26 

f n 







I 4 

1 

52 

Pa 






H 

Ul 
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6 


l ihe 

AOOR 

LAHL 

bl 

02 

HJ 

6a FfcHOH 

SOURCE 



1 17 

0 7i)C 


A3 

FS 

CD 

cr 

DATA 

Ht43tF5iCl>* 

CC.Flt3E*B3.4F, 13,00* 




FI 

3E 

H3 

4F 







13 

no 






11H 

0716 

0716 

01 

0? 

04 

OA 

BCD DATA 

Hi 01 1 02 i 04 ? 

OQ * 16 »32»64 , 00 * 




16 


64 

00 




It* 

071E 

071E 

C13 




SINE GAIA 

H* 03* 


1 20 








121 







GRG 

H i 75Q 1 


122 








123 







« INIALIZATIOn 

OF THE instrument • 

124 







* PARAMETERS AFFECTED incluoej * 

1 2b 







** (1) PR. status lOWER/UPPER * 

126 







* (2) CONTENT 9 

OF REGS* C4D 

ON CPU BOARD * 

127 








128 







44oo INITIAUZF 

o l«« 

1 29 

U7S0 


76 

2 1 



PPSlt 

il 

INHIBIT INTERCEPTS 

UQ 

0/b2 


7b 

f E 



CPSL 

m ff 1 

CLEAR ALU 

131 

0?S4 


77 

0? 



PPSL 

COM 

LOGICAL COMPARISONS 

132 

0/S6 


04 

FF 



LOO I i RO 

NSAM 


1 33 

07 S8 


FO 




WRTU « RQ 


NO sample/hlo or MUX 

134 

07S9 


74 

FF 



CPbL 

h»FF ' 

CLEAR ALL 

US 







«« 



us 







»»»* START MaIm program 

« » a a 

U7 

uisa 

07SB 

IF 

07 

riO 


main H$rA,ON 

SAMP 

SAMPLE ALL Data CHANNELS 

u« 

07SE 


IF 

OA 

E6 


bsta,un 

ROTX 

correct hx Data 

139 

0?61 


IF 

0 A 

71 


HSTA,UN 

RoTy 

CORRECT HY Data 

uo 

0?6A 


JF 

OA 

09 


bsl a*um 

rtvEC 

compute mori?. vector 

1M 

0 76 7 


3F 

07 

/ 0 


HSTa.uN 

RICH 

CALC. HEADING 4 0/P 

U2 

0 7b A 


3F 

08 

97 


8StA,UH 

OuTa 

a/p roll ^pitch 

1 **3 

0760 


IF 

07 

bB 


8CTA ,UN 

MAIN 

LOOP FOREVER 

144 







4#ft04#4#0*0**0040#000404#<»tH»a<i40O«040««M*»i»0tt0U40'»* 

US 










| *6 







a 



u? 

0770 

07 70 

77 

0? 



*ICH PP$L 

Lob 


ua 

U7/2 


OC 

04 

/(] 


LODa.hO 

Hy2 


149 

077S 


€C 

04 

*E 


COMd ,R0 

HX2 


1 bO 

0776 


\o 

OH 

14 


0CTa*6T 

Cosy 


Ibl 

0770 


1 E 

fiH 

2H 


tiCTA f LT 

bjNY 


IS2 

0 7 7E 


nc 

04 

71 


LODa # RO 

Hy2.1 

L& BITS 

1 S3 

U 7h 1 


EC 

1)4 

6F 


COMA # RO 

hk2* 1 



07ti4 


10 

OH 

14 


BC T * | (jT 

CnSr 


US 

0787 


JF 

OH 

28 


UCT a ,UN 

SINY 


U6 










LO 

Ln 


BEPRODtrcMLmr or the 

ORIGINAL PAGE IS POOR 



line 

lbd 

159 

160 
1 6 i 
162 
16J 
164 
16b 
166 
167 
loti 
lo9 
170 

m 

m 

173 

1/4 

175 

1 76 

\n 
i?a 
l 79 
ldO 
1U1 
1 ti2 
163 
lti4 

US 

166 

167 

168 

189 

190 

191 

192 

193 

194 
1 9b 
196 
19/ 
I9fl 
199 
500 
201 
202 
203 
2U4 
/Ob 
206 
2U/ 
20b 
209 


HIP ASbFHMLFR VFRMON SCO LE^El \ HEADING INSTRUMENT aSS£m L Y PROGRAM 197b PAGE 
ADDR LA-.L 6 1 62 13 Ha rRROH SOURCE 

OkG H » 7rt0 1 

»* A bOtiWoUllNE TO SAMPLE DATA f ROM ALL SENIORS ° 
« RtGiSlE* "0" CONTROLS SAHPUt/HULO C MUX * 
« tins 0.1.L2 SELECT hX.HY .HZ, THF T**PHJ * 

» DATA IS STORED IN DOUBLE PRECISION FOR'1 * 
* BEGINNING aT "DATA 11 *ITh HXCMSB) FOLLOWED RY * 
•» nX<LbB) * HY * M7» P I TCH , POLL • * 
a (1) COfsVMHS A/D DATA TO SIGN maGNITUOF ° 
» (21 CHANGES SIGN OF Hf L Hi (OHlENTATIOn) * 
“ <31 CORRECTS tfik OFFSET ERROR * 
<* FXlT wlTH rORRFCTEB DATA IN TAtiLt "OATH 0 * 








OftUll 

BLblN SlIQrOuT I NF 

««»« 

0 /riO 

0 7 *t 0 

Ob 

Ff 



SAMP 

LODI ,rfl 

-1 

OaTA INDEX 

0 7n? 


06 

OH 




LODI iN2 

5 

LOOP COUNTER 

0?ti4 


0/ 

01 




L001 ,R3 

l 

1ST sample 

0706 


75 

0 9 




CPSL 

C*wC 


o7fta 

07<»0 

F 1 




MtXT 

vtU Tu,P3 


SELECT bAMPi £ DATA 

U7ti9 


67 

BO 




i on i .in 

Hi 8U 1 

READY to hold 

0 7OR 


F3 





WPTU.R3 


hold oata 

0 7riC 


04 

HO 




LODI .HO 

MiftO* 


0 7b£ 


“0 





WHTC.HO 


RESET A/D 

0 7tiF 


20 





E0R2 

«o 


O7L0 


MO 





wrK.RO 


S f ART A/0 

0 7C 1 

0 ?C 1 

70 




TEST 

redd.ro 


READ LSd 

Ole? 


F4 

0 1 




ThI.ro 

HlOl * 


U7CA 


l A 

?H 




HCTm.n 

Tf ST 

BRANCH BACK IF FOCsl 

(J ?C6 


?4 

ff 




eori.ro 

H»FF 1 


o?ca 


44 

FU 




AnD1.ro 

UpPr 

STRIP OFF 4 LS HITS 

0 ?CA 


CO 

24 

SI 



STRa.ro 

DA T a ♦ 1 . H 1 . ♦ 

S 1 ORE LS 1/J 

0 7 CO 


30 





REDC.kO 


READ Mb 2/3 

0 /CF 


24 

T F 




EUR I .HO 

mff ♦ 


0 7 IJO 


CL) 

24 

4F 



STRA.HO 

UaTA-I . p 1 ♦ ♦ 

STORE MSB 

U703 


£»/ 

0 1 




adoi ,m 

1 

POINT TO NEXT DATA 

0 7i)5 


a7 

OF 




AMO l |R3 

LnwR 

POINT TO CHAM „ ONLY 

o7i)7 


FA 

SF 




HDRR.R2 

NeXT 

Tt-ST LOOP COUNTER 















CONVERT TO 

bjON MAGNITUDE 

0 7l)9 


77 

OH 




PPSL 

wC 

ENAHLf CARRY 

0/08 


06 

FF 




LUD1.R2 

“l 

iNpEX 

O7of) 


07 

0 7 




LODI »P3 

7 

LOOP COUNTER 

0 7 0 F 

n/fjF 

oE 

24 

SO 

- 

CONT 

LODA.RO 

DaTa.r2»* 

GET MbU OF DATA 

07E2 


IA 

1 1 




UCTR.N 

NSlU 

UK. IF SIGN*| 

U7F4 


Cl 





S UU 

R 1 


O/Fb 


2u 





tmi 

kii 

CLEAR HU 

07t6 


77 

01 




PPSL 

c 

SLT UP FOR SIJRT . 

0 7£R 


aE 

?4 

SO 



SU6A.R0 

Ua TA »R2 ♦ ♦ 

SUB » LS BITS 

07EM 


C£ 

64 

SA 



STRa.ro 

UaTM.PS 

STORE NF W OATA 

0 7tE 


20 





eokz 

Rn 


07lF 


Al 





SDb/ 

R 1 

SUB, Mb HITS 

07F 0 


64 

40 




I0R1 iRO 

H» Bo » 

SET SJGN^l I MF Ci ) 


136 



PIP AW>WUU VERSION &rU l£Vtl 1 HEADING INSTALMENT aSSEMLY PROGRAM i S>76. PAGE a 


LINE 

AODR 

LAOL 

Hi 

6? 

HJ B4 £ RROR 

SOURCE 



?id 

07F2 


CE 

64 

59 

STRA ,RO 

OATM-I.R2 

STORE MS BITS 

El] 

07FS 


1ft 

OH 


BCTH.UN 

F IN I 


212 

07F7 

0?F 7 

44 

7F 


RSTO AND l t HO 

H i 7 F ■ 

SIGN on - 0 

213 

U?F9 


C£ 

64 

SA 

stra , wo 

DaTm * R2 

STORE NEW DATA 

214 

07FC 


oE 

24 

SO 

LDDA , RO 

Data »N2* * 

GET Lb DATA 

21b 

U7FF 


CE 

64 

5A 

STRA ,RO 

OaTM* R2 

STORE AS NEW DATA 

216 

uao 2 

nans 

fa 

5ft 


F INI 0DRR ,W3 

CoNT 


217 






» CHANGEBSIGN 

Of Hy 6H2 < SENSOR ORIENTATION) 

216 

0604 


<)C 

04 

5C 

LODA * HO 

OftTM-2 

change sign of hy & hi 

219 

oao7 


no 

04 

5E 

LODa.HI 

OaTM *4 


220 

uaoA 


24 

<Jf) 


LORI ,R0 



221 

oauc 


?b 

Kn 


EOR1 *Hl 

Hi QO • 


222 

fl«0E 


CC 

04 

SC 

STRA t MO 

DAT 4*? 


223 

0611 


CD 

04 

bE 

STH A *R1 

DaTm+4 


224 






* COHHLCf FOR 

offset error 

IN hX,ny. «. HZ channels 

223 

0614 


0/ 

00 


LODI ,RJ 

0 

DO HX 1ST 

226 

0616 


3f 

OH 

2A 

BSTa,UN 

ofst 


227 

06 19 


0 7 

(IP 


LODI *H3 

2 

HY 

22 a 

uam 


3F 

u8 

?A 

B$TA,UN 

ofst 


229 

UttlE 


07 

04 


LOOI.R3 

4 

H2 

2 JO 

<1620 


If 

HH 

? A 

BSTA.UN 

ofst 


2J1 

0823 


CO 



NOP 



2J2 

0624 


07 

07 


LOO] *P3 

7 


233 






* CORRECT FOP 

ORTHOGONAL 1 TY ERROR IN HA AXtS 

234 

0626 


3F 

OF 

ft 8 

USTa,UN 

ORTH 


235 

0829 


17 



RE TL ♦ UN 



2J6 








137 



PJH ASSFMALfH VfwSlON SCO LtVLI \ hFADjNG iNSTRtMf-NT AS^t^LY PROGRAM 1976 nAGF 
L i»Mt AlJDP I AHL til dp Hj iU 1-fakOri SOURCE 




2Ja 

2J9 

2AO 0 B2 A OB^A 
2a l utipn 
2A2 Od?E 
<***3 l)B3fl 
“2AA UB3? 

2sb oajs 

?A6 UB3H 
2A/ dftJH 
2Aft utjJF 
2A9 Ubnl 
250 0Ba4 
2b 1 UtiA7 
2b2 ua^A 
2b3 (ihAO 
2bA UtsbO 
2b5 0852 

2b6 OUbA 

257 06S7 

2b a 0bb9 

2b9 Ub^C 

260 ObbF 

261 0«62 
262 (JftftS 
26J 

26A 

26 b 

2t>6 

2o7 

26H 

269 

2/0 

271 

2/2 0U66 UBS6 

2/3 OShft 

2/A 0869 

2/b 0B6C 

276 086F 

2// 08/2 

2/B 08/A U67A 

279 0677 ob/7 

2d0 

2«M Ob 7ft 

do? Ob/9 

2oJ 06/A 

26A Oa?C ub/C 

265 06 n fin ft 

tub, 0881 

2*7 OBdA 

268 Ob<i 7 

269 


OF 6A bA 

U 

AA ff 
Ab tM) 

CC 11 A 77 
ro oa /f 

OF 6 » SH 
CC UA /II 
UF 6 A Ofc 
CC 04 6b 
<*f 6A ftf 
CC UA 66 
OF 4/* l iA 

cc u<* /n 

/3 0 2 

7/ 10 

|F 09 98 
75 JO 

nC oa n 

( f 64 bA 
UC OA /ft 
<f toA Sft 
17 




LOb««*JO 

ua rM,pj 

STP2 

pi 

Anui »K0 

HP /f • 

AND I ,H\ 

h* 60 » 

STKA.kfl 

PSLT 

S T RA t H I 

SPtL 

LOUw,M> 

OATMH ,03 

STPA.HO 

Pstro 

LOU w . R(j 

ose.t,.o 

STPA.R0 

temp 

LUOA.kO 

user* l *Ri 

SlPA,Po 

7fmp*i 

LOCA.no 

SnST 

STPa.PO 

STCH 

CPbL 

Com 

PPSL 

PS 

BSTa#UN 

Sadd 

cpsl 

PS 

LOUw.no 

psl r 

SIHA.DO 

DATH,P3 

LODw.PO 

psltm 

stpa.po 

Pfc.TC.UN 

UaTMO ,P3 


MS 8vr£ OF DATA 

bUlIP SU>*’ 

SAVf- SIGN 

LS UY1C 01- JArA 
Hb HYTt OF nfrSFT 


bank*i 

BAnk-0 


b? 

20 
DO 
hb r»D 
05 02 

uU 6A 7 o 
ftD /) A /A 

< 0 <«. 7t, 

F 9 7b 


" r ™ ~ — 

H I M-uRY 'JlT. FOk A TWO-bYfE I NTEGtk 
MULI. flr A SlNG| E BYTE I NTCuLH 
HUL1 IPUIC9 JS IM OPWl ♦ 1 
MULIJMLTCaNQ IS IN DPR2.0PP2M 
DCSlLT wILL OF IN PSLf *BSIT*] 

(2 aylFS APF DISCARDED } 







«»o# 

df GIN Si JHnOii l INI 

9 ** « 

Oft 


bMPY 

PPSL 

wr 

SET MUDfc 




tim/ 

wn 


OA 

n 


STPA.PO 

PSLT 

CLEAP RESULT 

OA 

/a 


STPA.PO 

PSLIM 

CLEAH MtSuL r ♦ i 

OA 

79 


STp<a,hO 

(JPR I 

Oft 

/ A 


LODI ,K3 

a 

IUAD COUNT 

OA 

LOCO 

LODa , k2 

DPR t + t 

CtT MULTI PL I FL 



LUOO 

PHR.R2 


Rdt » hIOht a \ th carry 



*• THIb WOTATCS 

MDi. T IPL IFft 

tiY l HIT 10 GFT L5B INTU 




FOP^ 

PU 

CLEAP PO 


LUC l 
L0C2 


6 PL • HO 
nOPk ,kO 
LODI ,»U 
LOUA.RO 
ADDA, DO 
‘HPA.kO 
HO UP, PI 


LDCA 

2 

PCLI-1 tPl 
DPH2-1 ,«1 

PSL r- I * H 1 
L0C2 


GLT CARRY INTI* l SH 
fiNANCh* If C = 0 
INDEX 

ADO MULTIPLICAND TO PRODUCT 


FlMIS-i Inf Afin 


i 

I 


£ 

> 

§ 

j-w 

CO 

O 

o 

w 


CO 

CO 


JS* 1 ® OP THE 



RIH ASSFM-LFR VeKSl »N S fU Lf Vtl l HEADING iNSTRt 


•Mr-Nl A5^b M Ly pPOGRAH 1976 RAr>£ lu 


L IfNt 


AU|)R l A \L ul *V iO ti c* (-W^UU SOURCE 


2V0 

0069 

0bR9 

Ob 

FE 


291 

i) 68 « 

OHHB 

no 

67 

79 

ay 2 

OUbE 


so 



293 

liHtiF 


CO 

61 

79 

294 

u892 


Ij9 

7 7 


?95 

UB94 


Fa 

M 


296 






297 

0B96 


W 



2va 






?99 






300 






301 






302 






303 

0897 

na9 7 

oc 

06 

nO 

706 

U 89 A 


C3 



305 

U&9H 


6 4 

nn 


JG6 

U69U 


1 A 

14 


30? 

0B9F 


2 C 

0 4 

62 

30« 

0dA2 


1 A 

06 


309 

116 AA 


06 

07 


710 

Ob A A 


04 

in 


ill 

OOAO 


IB 

16 


312 

UbAA 

ftbAA 

?0 



313 

06Afl 


04 

10 


714 

OBAO 


IB 

UT 


315 

obaf 

A8Af 

?C 

04 

aa 

116 

08b? 


19 

06 


31 / 

UbfiA 


ft 4 

ft i 


31 a 

0 brt 6 


U4 

in 


319 

OHim 


Id 

ft ♦ 


320 

UbfcJ A 

08HA 

04 

1 ) rt 


321 

OBijC 


04 

10 


322 

06bf. 

nH^E 

7F 

OR 

CC 

323 

06C1 


ns 

04 


326 

0hC7 


OF 

04 

62 

325 

06C6 


IF 

04 

CC 

326 

08C9 


US 

OH 


327 

oocn 


l t 



326 






729 

08cr 

rmcc 

4 / 

n 


JJU 

OdCb 


CF 

34 

/C 

3Ji 

0 6UL 


?0 



332 

0 Wl )2 


cc 

0 4 

7B 

337 

0bu5 


ns 

s.\ 


336 

Ub|)/ 


CC 

in 

/ A 

335 

ObuA 


if 

im 

46 

3 36 

Ob 00 


tvc 

0 4 

/a 

337 

OdFO 


JF 

OC 

44 

3JR 

uatT 


?S 

Ff 


.139 

oats 


17 



160 







LOC4 

L0C5 




l_OOl.Nl -? ROT. THE PRODUCT 

LODA.RO RSLT~ 2 S 4 * 2 *Ri 

ftNR.AO 

STftA.WO RSLT-2S6+2.R1 
BXRR.Rl L(\C5 

8DRR* R3 LoOft FINISH THE LOOP 

♦ **»***•»* ***»*'**»«*»♦<»» 


RETC.UN ^ ^ 


* A SUdPQuTl <€1 Tc\ OUTPUT PITCH L ROLL OAtA 

* CALLtO BY MAIN PROG • , FETCHES CONVERTED DATA* 

* FROM TAQLE j, OAtM‘' then o/p*s acu angle * SIGN* 


OUTA 

LOOA.RO 

OaT A*f> 


STR 2 

K7 


AnOI.RO 

h i a o * 


8 CTR.N 

NTH 


EOHA^NO 

Oa Tm *8 


BCTR.’N 

NPH1 


LOOl.UO 

7 


tfftlt.NO 

SRI T 


(JCtH.UN 

AGL 

NPH 1 

EUK 2 

K A 


WRTh.RO 

SPIT 


BC7H.UN 

AGL 

NTH 

EORa.ro 

OaTH +8 


&cn<,p 

NPH 


LOOl.PO 

3 


w«rt,wo 

SrUT 


BOTH t UN 

AGL 

NPH 

LOU i ,ftG 

a 


WRIt.UQ 

SPUT 

AUL 

bsta.un 

LOU 


WRTE.Rl 

1h T a 


LOU* .1<3 

Da T M * a 


BS T A* UN 

LOD 


wRTt.Rl 

RFTt.UN 

Pm! 



LOD 

AND J »R3 

H.7F* 


STH* »f<3 

UpR 2 # 1 


Lu Hi 

R^ 


STRa.RO 

c pr 2 


LODI.NO 

Hi Sa * 


SIRa.RO 

rpRI H 


HSTA.UN 

SmPY 


LODA.RU 

RSLT * 1 


flSTA.UN 

BCOA 


LOR I ♦ R 1 
Rf TC.UN 

HiFF • 




FtTCH THETA ImSRI 

save SIGN 

BR* ON MEG THFT A 

COMPARE signs 


BOTH POSITIVE 


PHI NFG. 


PITCH NtG.6 ROLL POS. 


BOTH NEGATIVE 
F UR# ANGLE (GcD) 
OUTPUT PI TCH 
FETCH roll 
angle <bcu> 

OUTPUT PHI 


MULTIPLICAND 

CLEAR LS rfYTF 
F ACTOr 1 . 35 1 S6) 
MULT I RL IER 
BINARY AnGLF 

angle 

ANGLE 4 toCD > 




H 

Oj 

LO 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR, 



unit 

36 1 
3** <* 
16b 
iHh 
)W 
16R 
3^9 
Isu 

lb l 
3b? 
IbJ 
3S4 
3bS 
Jb6 
3b7 
JbH 
3b9 
160 
)6l 
3b<> 
163 

36 *♦ 

165 
366 
3b / 
i6B 
Ib9 
3/0 

JM 

j72 

in 

1/b 
3 /<> 

)/7 

37 a 

179 

PjO 

Jdl 

362 

363 
JdA 
36b 
366 
16/ 
360 
369 

390 

391 
Jv? 
19J 


HIP A6SPMMLfP VfPSjON SfU LFVtl 1 hfAUlNG I t^STPf iMFNT ASSFMLY PROGRAM 1976 PAf,t 11 
AUGR LA fc <L til n? i2 bA F PPOK SOURCE 

» S*K. kCTX TO CoWWECT Hx'OAlA USING ALGORITHM “ 

* MAsnX* *C0S t* a ITCH) ♦ S I N <P l fCHJ MHYm«SIN( ° 


<* put Li *«ZM‘>r0S(H0LL> ) * 

» tWftP A|lh “Paw TATA" IN ThBLE UATM * 

“ txll wlTri *0 )IF ItO HX fN "uATX" * 

» CALLS S *H • H7HI * 

* CALLS S.H.S, »SAO0“ k *'SMPt" TO 00 « 

* SlGNtU MAGMtUoE AOO/bllBT a NO MULTIPLY * 




06L6 

Ofirb 7? 

<ib 


|<0TX PHbL 

wr 

Af< { TH WITH CAUbY « 

UdEP 

oC 

OA 

ftO 

LOUA,H0 

0 A t M * ft 

The r A (PITCH) 

OBEB 

ri 



STkZ 

^1 


OBEC 

r i 



STPZ 

HI 


OBtl) 

A 5 

n 


* AN 13 1 *P 1 

Hi 7f • 

S/M1P OFF SI69 

ObLF 

A/ 

bn 


A\iU I »w3 

Hi 60* 

SAVh SIGN 

ObFl 

/ f 

OA 

ft 7 

strati 

Si NO 

SIGN T ML T A 

06F4 

oti 

6S 

00 

LOU A ,^0 

srotRi 

S 1 N { ThE T A ) 

OOF 7 

( c 

OA 

66 

b T p A t K 0 

SjNDM 


OOF A 

t*V 

6H 

69 

LOUM «ND 

cns^fti 

CUS( TntT Al 

ObFO 

Cl 

OA 

ftS 

Sfpu^o 

CnSf 


09(10 

( c 

OA 

7A 

STM* ,«0 

OPHI * \ 

mult ihliep 

0903 

ii c 

OA 

S A 

LUU*» .*0 

UaTm 

HA Mbrtl TS 

0906 

Cl 



S T H/ 

HI 


090 7 

aa 

1\ 


AWUI *H 0 

H 1 7 f* 

STRIP OFF S l Go 

0909 

*.b 

60 


ANU1 **] 

h.diH 

save sign 

09 oR 

(f 

it A 

7d 

srf<A,po 

OP US 

mult IPlICANO 

090E 

f J 

0 A 

Mb 

SlHA»9l 

STOa 


09 1 L 

JC 

0 *♦ 

sb 

LuiJA *P0 

b 6 1 * * 1 


09 l A 

( c 

OA 

/C 

STH*,PO 

0PR?M 


0917 

if 

00 

66 

BS TA *UN 

SmPy 

f Ohm h x “ COS (THE 

09 1 A 

OC 

OA 

n 

LOU A *110 

KSLT 

MOVE HPOOuCr 

09\f) 

cc 

ft A 

H|l 

ST«A f HO 

cruf 


09?(1 

nc 

OA 

7/1 

LOU-.bO 

«SLfM 


09^3 

LC 

ft A 

HC 

bT9A f W0 

O.iUF 


U9?6 

jf 

Of 

aO 

bst 0 # un 



09^9 

iiC 

OA 

ft a 

LUUA ,kO 

S | N J ♦ l 

SIN (THETA) 

09^C 

t C 

OA 

/A 

sr«A f uo 

OPH 1 ♦ 1 

mult iplier 

092F 

If 

Oft 

66 

dST A f UN 

S» HY 

hZ* 1 ^ In ( ThETaJ 

093? 

»( 

OA 

6 7 

LOO* ,Po 

b l N'i 

SIGN OF TmETa 

09 JS 

CC 

OA 

/ t. 

STOA.t'U 

ScEb 

Fok Us£ By SaOU 

093H 

o 0 

OA 

a 6 

LOUA ,rt 1 

SlGx 


u9 JH 

( u 

OA 

70 

b THa , PI 

StEh 


09 IE 

hC 

ft A 

Hb 

Li)OA # R0 

CiiUf 


09A 1 

ou 

OA 

nc 

LOO A ,H J 

OMUI- 


09AA 

U 

t'A 

6b 

Si H A 4 P 0 

1FMP 


U9a 7 

ru 

OA 

*>6 

STPA,P1 

TFHM* 1 


09 a A 

/b 

ft? 


CPSl 

COM 

FLAG a AOO 

09aC 

tf 

09 

96 

HST \ ,l«N 

SaUO 

FUpm wtw HX 

09Af 

OC 

0 A 

76 

LOOA *po 

HSLl ♦ 1 

HOVE Ntw riX 

U9b? 

( c 

ft A 

60 

sroA f pu 

f/A 1 A * 1 


u9bS 

OC 

OA 

77 

LOOA.NQ 

halt 



M 

O 


®S§BUGffiILrFY OF THE 
OB W& MSB IS POOR 



1976 PAfiF 12 


LI Alt 

394 

396 

3 96 

397 
390 

399 

400 

401 
41) 2 
41)3 
4U4 
4Ub 

406 

407 
4ua 

409 

410 

411 

412 

413 

4 14 
4 16 
416 

4 1/ 
410 
4 19 

420 

421 

422 

423 

424 

425 
w26 
42/ 
423 
42 / 


^IP A^S^HLPR VFQSI'IN ScU LFVti 1 HEADING livSTR'lMPNT 3SSJ-HLY PROGRAM 
A(Jt)W L AHL 111 B2 ri 1 tiu FRRUR SOURCE 

DaTX 


09btt 

OvSH 


CC 04 ff 
17 


04 bC 

096 E 

0 V 6 F 

09 O 2 

0965 

0967 

0969 
U96C 
096F 

0970 

0971 
U972 
097 1 
0975 
09/a 
09 7 a 
09 /£ 
U931 
U9ri4 
096/ 
09rfA 
U9IH 
09 hF 
09 y l 
0992 
u9>6 

099/ 


09SC 77 
?0 
LC 
CC 
07 
4b 
CO 
CE 

096F 51 

62 
20 
no 
FB 
oc 
«c 

CC 

oc 
ac 
1 F 
OC 
50 
CC 
nC 
50 
CC 
M 

17 


09m7 


oHO»o»n(Hfitaui)ita» 


STRA . wO 

RETC.UN 

(> O O <* » « U O O (> i> 3 u o t* » iJ _ 

•* suukOuT inl ro squaw two ^ byte values 
* enter with value in pkmsb> 4 h2<lsb) « 

“ EXIT WITH 2 MSd In WSLT RSlT+ 1 a 




o««» 

BEGIN 



0 « 


sou 

PPSL 

wr 





tori 2 

Rn 

“ * 

04 

77 


STHa.ro 

HSLT 

CLEAR TEMP STORAGE 

04 

7 rt 


STRa.ro 

H 5 LTM 


10 



LOO I *k 3 

16 

COUNTER 

7 F 



an!) I ,ki 

H » 7 F 1 


04 

/ri 


STRa.ri 

OPR? 


04 

7 C 


STRA.R 2 

OPR 2 *| 




ROTA 

RRR . ft I 


ROTATE MULTIPLIER 




RHH.R 2 






E 0 R 2 

Ro 





RRL *R 0 


GET LS dll INTO RO 

IP 



BOHR, NO 

noao 

NO A 01 ) IF C »0 

04 

/fa 


LODa.HO 

rslt*) 


04 

/C 


AOOa * ftO 

OpR 2+ 1 


04 

70 


S TftA , RO 

Rstr-i 


04 

77 


LUUA.kO 

RsL r 


0 4 

7)1 


ADOA.RO 

UPR? 


09 

A A 


BCTrftUN 

NO AO* 3 


04 

n 

NO A 0 

LOUA.KO 

RSLT 




** 

HRR.RO 



0 4 

// 


STRA.hO 

HSLT 


1)4 

/ii 


LOOa.ro 

RSLTM 





RRR .rin 


RUTATh LS 1/2 KESUI 

< it 

/fa 


SIha.ro 

RsLT ♦ 1 


SH 



H 0 RR.H 3 

Rn T ^ 

CONI „ MUL f * IF Of Q* 


RFTu.UN 


END Oh 6UHROUT INF 


© W 

S w 

g 2 
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APPENDIX C 


The transcendental functions used throughout the heading 
computation algorithm were implemented using a table look up 
procedure. To generate the respective look up tables in com- 
puter memory data was first generated using algol programs. 
This technique expedited modifications to tabular data and 
provided output data in a convenient (hexadecimal) format. 

Programs that calculated Cos (0) and Cos 2 (0)to eight 
bit and sixteen bit resolution respectively are included. 
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INTRODUCTION 


The fact that the Earth has a surrounding electric field has been 
known for centuries, but the study of the effect of clouds on this 
electric field is a fairly new area of study. Ideally the earth's 
electric field is a perfect vertical, however since clouds are not 
electrically neutral, their presence distorts the earth's field; 
nearly neutral clouds causing slight perturbations while heavily 
charged clouds actually reverse the direction of the resultant field. 
The most highly charged clouds belong to the cumulonimbus family. 

The report evaluates various electrical models of cumulonimbus clouds. 
The resultant field of each model is compared with actual readings of 
the earth's electric field in the presence of a cumulonimbus cloud 
at various heights and distances from the cloud. Using the actual 
field readings, this report will develop a new electrical model of a 
cumulonimbus cloud. 

This paper also reviews the electrical properties of the atmosphere, 
particularly of clouds and suggests some possible uses of the 
knowledge of the field around a cumulonimbus cloud. 
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II. ATMOSPHERIC ELECTRICITY 

Atmospheric electricity is defined by Dolezalek [1] as "the science of 
eveiything electric in the atmosphere between the solid and- liquid 
surface of the earth to the lower regions of the ionosphere, including 

r 

the boundary layers of both." 

Like every science, atmospheric electricity has its own sign conventions. 
Distance is measured positively upward. Current flow to the earth is 
.considered positive. Electric potential is measured relative to earth. 
The terms "field" and "potential gradient" are in general used to 
denote the same quantity, i.e. the partial derivative of voltage with 
respect to height. The units used are the MKSA system. 

1 All lines of force commence on a positive charge and end on a negative 
charge. The density of lines of force across any area gives a measure 
of the field. If the line of force terminates on the earth's surface, 

- it must enter vertically. Any change in space charge will immediately 
effect the lines of force and thus, the electric field. 

As space charge appears, a potential gradient is produced immediately 
at the ground. The air-earth current will change as will the lines 
of force. To determine the. potential gradient change produced due to 
the appearance of the space charge, one must know the relaxation time, 
e/X, of the air where e is the permittivity of the air and X is the 
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conductivity. Near the earth's surface the relaxation time generally 
ranges between 5 to 40 minutes depending to a great extent on the 
pollution level. The higher the pollution level, the lower the 
conauctivity and, therefore, the higher the relaxation time. At an 
altitude of 18km the relaxation time is about 4 seconds while at 70km 

I 

it is on the order of 10~® seconds. At the earth's surface it is 
10' 6 seconds or less. The potential gradient at some time after a 
change in space charge can be given by the equation [2] 

E = E-j e" t/T + E 2 (1 - e“ t/T ) (1) 

where is the potential gradient immediately after the change, E^ 
is the potential gradient when the conditions are ac.'in quasistatic, 
t is the relaxation time and E is the potential gradient in the 
intervening time. 

Ions are relatively rare when compared to the total "'..mber of molecules 
in the atmosphere. A cubic meter of clean land air *ill contain 

oc 

approximately 800 million ions. out of 10 molecules. The conductivity 
of the atmosphere may be written as [1] 

X - Zn, + e f + k* + Sn 7 e' kj" (2; 

where n. + is the number of positive ions with a charge of e* a-d a 
mobility of k + , n." is the number of negative ions * ; :h a chs 

* O 

el* and a mobility of kT , and X is the conductivity The mc.~' 

J J 



is computed from the potential gradient E and the velocity of the 
particle v by 

k = v / E (3) 

or the velocity acquired in a potential gradient of 1 volt/meter. 

i 

The mobility is highly dependent on physical properties of the ion 
such as size and shape. 

Since the air-earth current is continually bringing positive charges to 
the earth there must be some restoring component which will complete 
the circuit carrying the positive charges upward again. The most 
popular theory on the generation and continuation of the electrical 
characteristics of the atmosphere is the condenser theory. [1,3] 

In this theory the ionosphere is positively charged with respect to 
the earth. The air-earth current tends to continually discharge 
this global condenser. However, this process is compensated by 
thunderstorms which act as the generator restoring positive charge 
to the ionosphere. Figure 1 shows the schematic of this global circuit. 
The lower portion of a thunderstorm is negatively charged while the 
upper portion is positive. The world wide occurrence of thunderstorms 
also corresponds to the diurnal variation of the earth's electric field. 

The study of atmospheric electricity is divided into two quite different 
divisions, fair weather and disturbed weather. Disturbed weather is 
defined by Dolezak as "when we have hydrometers in the atmosphere - 
fog, rain, snow - or when we have high winds, in particular in 
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connection with blown-up dust or snow, or when we have much of the sky 

covered by substantial clouds." [1] The absence of these phenomena is 

considered- fair weather. Fair weather conditions are predictable and 

one can assume a quasistatic state to determine its properties, while 

disturbed weather condi tons can alter -rapidly and a quasistatic state 

< 

cannot be applied. 

Fair weather conditions have several characteristic quantities. Voltage 

increases with height. The potential gradient is positive and essentially 

vertical. It is constant for the first few meters and shows a progressive 

decrease with height at approximately 100 meters. The conductivities 

are also constant for the first few meters, then show a marked increase 

with height at approximately 100 mei: o rs due to an increase with altitude 

of ionization by cosmic rays and a decrease in pollution. The current 

density is the same at all levels and is equal to the voltage of the 

2 

ionosphere with respect to earth divided by the resistance of a lm 
column of air from the earth to the ionosphere. Tables 1 and 2 list 
fair weather parameters as compiled by Dolezelek. [1] 

The disturbed weather phenomena is much harder to typify. The potential 
gradient can have horizontal variations of large magnitude. There are 
added currents carried by precipitation and lightning. Constant 
changes in conductivity occur due to fresh charges. It is in this 
still unsolved phenomena that much of the current work is being 
performed . 



Part of Atmos- 
phere for which 
the Values are 
Calculated 

Currents, I, in 
A; and current 
densities, i, 
in A/m 2 

Potential 
Differences, U, 
in V; field 
strength E in 
V/m 

Resistances, R, 
in ft; Columnar 
res . , R c , in ftm 2 
Resisti vi ti es , 
p , i n ftm 

Volume element 
at about sea 
level , one 
cubic meter 

i=3x!0" 12 

E =1 .2xl0 2 

P o =4xl0 13 

Lower column of 
1 nr cross 
section - sea- 
level to 2 km 

same as 
above 

at upper end: 
U -j =1 .8x1 0 5 

R cl =6xl0 16 

Volume element 
at about 2 km 
height, 1 m 3 

same as 
above 

e 2 =6.6xio ] 

P 2 =2.2x10 13 

Center column 
of 1 m^ cross 
section - 2 
to 12 km 

same as 
above 

at upper end: 

U =3.15x105 . 
m 

R =4.5x1 0 1 6 
cm 

Volume element 
at about 
12 km height, 

1 m 3 

same as 
above 

E ]2 =4.2x10° 

P 12 =1.3xl0 12 

Upper column 
of 1 m^ cross 
section - 12 
to 65 km 

same as 
above 

at upper end: 
U u =3.5xl0 5 

R_=l .5xl0^ 6 
cu 

Whole column 
of 1 nr cross 
section - 0 
to 65 km 

same as 
above 

at upper end: 
U =3.5x105 

R c =1.2xl0 17 

Total spher- 
ical capacitor 
area: 

5x1 0 1 4 m 2 

1=1 . 5xl0 3 

U = 3. 5xl0 5 

R =2.4xl0 2 


Table 1 
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Part of Atmos- 
phere for which 
the Values are 
Calculated 

Conductances, 6, 
in ; Col . 

conductances 
G , in ft“}m -2 ; 
total cond. 

A, in ft - }m - } 

Capacitances, C, 
in F; Col . 
capaci tances 
C c , in Fm -2 ; 
Capacitivities 
e, in Fm - } 

Relaxation Times 
t, in seconds 

Volume element 
at about sea 
level , one 
cubic meter 

A o =2.5xl0“ 14 

1 

e o =8.9xl0" 12 

t o =3.6x10 2 

Lower column of 
1 m 2 cross 
section - sea 
level to 2 km 

G cl =1.7xl0' 17 

C c1 =4.4xl0“ 15 

t c1 =2.6x10 2 

Volume element 
at about 2 km 
height, 1 m3 

A 2 =4.5x10" 14 

e 2 s 8.9x10“ 12 

t 2 =2x!0 2 

Center column 
of 1 m 2 cross 
section - 2 
to 12 km 

6 =5xlG~ 17 

cm 

C =8.8xl0“ 16 
cm 

t =1.8x1 O 1 
cm 

Volume element 
at about 12 km 
height, 

1 m3 

A-i 2 =4. 0x10“ 13 

e 12 =8.9xl0“ 12 

t-| 2 =1 .2x10^ 

Upper column 
of 1 nr cross 
section - 12 
to 65 km 

G =2.5xl0~ 17 
cu 

C cq -1.67x10- 16 

t cu =6.7xl0° 

Whole column 
of 1 m 2 cross 
section - 0 
to 65 km 

G =8.3x1 0" 1 8 

C c =1.36xl0“ 16 

t=1.64x10 1 

V 

Total spher- 
ical capacitor 
area: 

5xlQ} 4 m 2 

G=4.2xl0“ 3 

C=6.8xl0“ 2 

t =1 .64x1 0 1 


Table 2 
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III. METEOROLOGY 

The nomenclature for the various sections of the atmosphere are not 
agreed upon by all who study the atmosphere, however the divisions 
and their names as used by Dobson [4] are shown in Figure 2. 

i 

The lowest cloud in the troposphere is known as fog. This occurs 
when the ground is within a few degrees of the dew point. Fog can 
be produced three ways. Radiation fog is created on clear nights 
if the earth cools greatly since without a cloud cover the heat 
radiates into the upper atmosphere. A slight breeze will then bring 
the earth to the dew point. Advection fog is caused by warm moist 
air blowing over a colder surface and becoming chilled to its dew 
point. Frontal fog may be produced when cold air mixes with warm 
moist air. Other cloud types are summarized in Table 3. 

The cloud names come from four basic words: cirrus, cumulus, stratus 

and nimbus. Cirrus clouds are composed of ice crystals and are 
delicate and curly in appearance. Sun shines through them without a 
shadow. Cumulus clouds are lumpy or billowing forms. The sheetlike 
cloud layers are stratus clouds. Nimbus clouds are stratus or 
cumulus that develop a "head" or thickness. "Alto" is also used 
as a prefix to indicate clouds a.t intermediate heights. 

When two air masses of differing temperatures, pressures and relative 
humidities mix, a front is formed. The colder air mass sits 
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Figure 2 



n 


CLOUD NAME 


AVERAGE HEIGHT AND RANGE IN HEIGHT 


Cirrus 

Cirro-cumulus 

Cirro-stratus 


Alto-cumulus 

Alto-stratus 


Strato-cumulus 

Stratus 

Nimbo-stratus 


Cumulus 

Cumulo-nimbus 


High Clouds 


30.000 ft. 

20.000 ft. 

35.000 ft. 


(10,000 to 50,000 ft.) 
( 8,000 to 36,000 ft.) 
(15,000 to 45,000 ft.) 


Middle Clouds 


12.000 ft. ( 3,000 to 27,000 ft.) 

15.000 ft. ( 5,000 to 35,000 ft.) 


Low Clouds 


5.000 ft. 

2.000 ft. 
2,500 ft. 


( 1,000 to 15,000 ft.) 
( 50 to 6,000 ft.) 
( 200 to 18,000 ft.) 


Clouds with Vertical Development 

2,500 ft. ( 1 ,000 to 10,000 ft.) 
Tops may extend to 20,000 ft. 


2,500 ft. ( 500 to 10,000 ft.) 

Tops may extend to over 35,000 ft. 


Table 3 

(from Aerology for Pilots [5]) 
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wedgelike under the warmer air mass. The air that is displacing 
another air mass gives its name to the front; thus, if cold air is 
displacing warm air it is known as a cold front. Likewise, in a 
warm front warm air displaces the colder air. An occluded front 
exists when two colder air masses trap a third air mass between them 
forcing the third air mass aloft until it dissipates. Figure 3 shows 
cross sections of the warm and cold types of fronts and the clouds 
associated with them. 
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Altitude (TO 3 feet) 



Altitude (10 3 feet) 



COLD FRONT 
Figure 3 

(from Introduction to Aviation [6]) 
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IV. EXISTING MODELS 

The cumulonimbus cloud differs from other types of clouds in its rate 
of accumulation of charge. A collection of charge's, such as a cloud, 
in a conductive medium will attract an equal but oppositely charged 

i 

screening layer and, therefore, cause no appreciable change in the 
electric field. However, since the rate of increase of charge in a 
cumulonimbus cloud within the time constant of the surrounding- 
atmosphere approximately doubles the amount of charge [7], an 
effective screening layer can not be generated and, therefore, the 
cloud effects the neighboring electric field. Models of cumulonimbus 
clouds can eliminate the screening layer and simply use the net 
increase in charge as the only charge present, 

Experimantors have been attempting for years to gather data on the 
electric field in the vicinity of a thundercloud, i.e. a cumulonimbus 
cloud, in order to generate a valid electrical model for the cloud. 

The measurements were originally made at ground level but more recently 
include airborne measurements in, around and above thunderclouds. 

In 1948 O.H. Gish and G.R. Wait [8], in an effort to show that thunder- 
clouds could supply the negative current flow' necessary to maintain 
the general electrification of the earth, gathered measurements of 
the electric field strength -above thunderstorms. Their findings were 
as follows: Of 87 traverses across thunderstorms, in 22 data was 
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incomplete or unsatisfactory, 31 profiles corresponded well to a 
bipolar model while the remaining 34 were of a more complicated 
nature where a bipolar model was only a good first approximation, 
additional dipoles would be required for an accurate model. 

A model was proposed using the data gathered on a flight an an altitude 
of 43,000 feet on October 28, 1948. This data was typical of a storm 
when little lightning activity was present. The data and associated 
model are shown in Figure 4. A complete mapping of the electric 
field in an area 21km in horizontal distance from the cloud's charge 
concentrations and to a height of 16km is shown in Figure 5. 

Israel [9] lists two other proposed electrical models for a thundercloud, 
one from Simpson and Robinson and the other from D.J. Malan. The model 
proposed by Simpson and Robinson is often referred to as the "Classical 
Model." It consists of a positive charge of 24 coulombs at a height of 
6km, a charge of -20 coulombs at 3km and a charge of 4 coulombs at 
1.5 km. The Malan model has a charge of 40 coulombs at approximately 
10km, -40 coulombs at about 5km and 10 coulombs at about 2km. The 
electric field map for the Simpson and Robinson model is shown in 
Figure 6 and for the Malan model in Figure 7. 

The program used to generate the data used in Figures 5, 6 and 7 and 
the data listings from the program are contained in Appendix A. The 
charge system for each of the three models used as inputs to the 
program are shown in Figure 10, utilizing the method of electrical 
images as described in the next section. 
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POTENTIAL GRADIENT (VERTICAL COMPONENT) OVER THUNDERSTORM, 
OCTOBER 28, 1948, STORM NO. 1, TRAVERSE NO. 8, (A) OBSERVED, 
(B) CALCULATED FOR MODEL SHOWN 


Figure 4 

(from Gish and Wait [8]) 




A = 39 coulombs 

B = -39 coulombs 


ELECTRIC FIELD MAP FOR 
GISH AND WAIT MODEL 


Figure 5 




VERTICAL HEIGHT (km) 



HORIZONTAL DISTANCE (km) 


A = 24 coulombs 
B =-20 coulombs 
C = 4 coulombs 


ELECTRIC FIELD MAP FOR 


SIMPSON AND ROBINSON MODEL 
Figure 6 


HORIZONTAL DISTANCE (km) 


A = 40 coulombs 
B =- 4 0 coulombs 
C = 10 coulombs 


ELECTRIC FIELD MAP FOR 
MALAN MODEL 


Figure 7 
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V. GENERAL NEAR-FIELD MODEL 


The field at any point in space due to a single point charge Q can 
be computed from Guass's law, 

§? * ds = / v ? dV (4) 

where D is the electric flux density, V is the volume enclosed by 
the surface S and z. is the charge density within V. If one considers 
S to be a spherical surface with its center at Q (Q being the only 
charge wi + hin S) , the equation becomes 

ft = (Q / 4irr 2 ) ? (5) 

Since & = sE where e is the permittivity of the substance then 

E = (Q / 4itr 2 e) r (6) 

The direction of E is radially outward if Q is positive and radially 
inward if Q is negative. See Figure 8 where Q is a positive charge. 
The magnitude of E is inversely proportional to the square of the 
distance to Q. If the permittivity is constant over all space 
concerned, the field created at a point P due to more than a single 
point charge may be computed by the vector addition 


t = 


4ne 


(^J-2 r lp + ~-2 r 2p + 

r lP IP r 2P dy 


NP 


-2 r Np ) 


(7) 


where Qj is the charge located at a point J and rj p is the distance 
between Qj and P. The permittivity used for this problem is that of 
free space, i.e. 8.854 x 10~^ 2 farads per meter. 
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The earth, here considered flat over the area under consideration is 

an equi potential . To produce an equipotential surface one can use the 

method of electrical images, where an imaginary distribution of 

charges inside a conducting body of an electrostatic system is 
1 

determined which would produce exactly the same field outside the 
body as that produced by the induced- free charges over its surface. 

For the system of a conducting plane, one can look at a simplified 
problem of a single positive charge located above the plane. See 
Figure 9. As shown in the figure, an equal but oppositely charged 
point is placed an equal distance below the plane as the positive 
charge is above the plane. Now if the plane is removed, the field 
above the plane would remain unchanged. Therefore, the field produced 
by all the charges induced on the conducting plane may be reduced to 
the field produced by the single point, -Q. This method can be 
expanded to include an oppositely charged point below the plane for 
each charged point above the plane in a more complex system. The 
systems suggested as thundercloud electrical models by Simpson and 
Robinson, Gish and Wait, and Malan are shown in Figure 10 using the 
method of electrical images. 

Each of these models was derived using the data collected by the 
particular experimentor under a particular set of conditions. The 
models fit the individual data, yet they vary considerably in the 
overall field maps as shown in Figures 5, 6 and 7. For example. 
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Figure 9 
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at 500 meters horizontally from the charge centers and at gound level, 
the Simpson and Robinson model yields a field with a magnitude of 938 
volt/meter, the Gish and Wait model yields a magnitude of 26122 
volts/meter, while Mai an' s model has a magnitude of 19991 volts/meter 
at that location. I.f a model could be generated to best fit a 

r 

combination of various field readings taken at different heights 
and different times under a variety of conditions, the result would 
be a more universal model of a thundercloud. 

The model developed in this paper is found by starting with the 
charge center values and locations determined by Simpson and Robinson, 
Gish and Wait, and Malan. The field strength is calculated for 
selected locations using one model at a time and compared to actual 
values measured by various experimentors . The measured values used 
for the comparison and the source of each value are shown in 
Table 4. The readings from Wormell [10] are average potential 
gradients taken over all directions from the charge centers at 
ground level immediately before a lightning discharge. 

The calculated value of the field using Equation 7 at location I, 
EFLD(I), is subtracted from the measured value, EF(I), and the 
percentage of variation computed as 


EF(I) - EFLD(I) 

efTT) 


100 


( 8 ) 



26 


DISTANCE FROM 
CHARGE CENTERS 
(METERS) 

MEASURED 
FIELD STRENGTH* ** 
(VOLTS/M) 

SOURCE 

1 2500 

- 1310 

Wormell [10] 

2 5000 

- 1160 

It 

3 7500 

- 760 

It 

4 10000 

- 460 

II 

5 12500 

- 260 

It 

6 15000 

- 135 

1) 

7 17500 

- 60 

II 

8 20000 

0 

II 

9 13100 

33000 

Gish & Wait [ 8] 


Note: 1 through 3 are horizontal distance with measurements made at 

ground level. 9 is a vertical height directly above the charge 
centers . 


*due to cloud only. 


Table 4 
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The resulting value is then squared. The square of the percentage of 
variation is -computed for each location in Table 4 and added together 
to produce -a quantity ESQRT. The amount of charge at one of the charge 
centers is now varied in increments of .01 coulombs until a minimum is 
found for ESQRT. The process of calculating the fields and computing 

f 

ESQRT is begun again with this new charge quantity replacing the old 
value for the selected charge center and a different charge center 
quantity is varied until a new minimum is found for ESQRT. The 
program continues to vary one charge center quantity at a time, each 
time finding the minimum ESQRT, until ESQRT is less than or equal to 
100 or until ESQRT reaches a local minimum. The program listing and 
resulting output from the program are contained in Appendix B. 

The resulting field maps of the Simpson and Robinson; Gish and Wait, 

and Malan models, after being modified by the program shown in 

Appendix B, are shown in Figures 11, 12 and 13 respectively. The 

progressive evaluation of the charge center values are shown in the 

output listings in Appendix B. The charge center values used for 

the modified Simpson and Robinson model are 4.6 coulombs at 6km, 

-11.3 coulombs at 3km and 5.97 coulombs at 1.5km. ESQRT had a value 
2 

of 5.22 x 10 at these values. The charge centers for the modified 
Gish and Wait model are 1.52 coulombs at 6096 meters and -2.64 coulombs 

O 

at 3048 meters; ESQRT equals 6.24 x 10 . The charge centers for the 
modified Malan model are 20.1 coulombs, -36.8 coulombs and 12.7 coulombs 
at 10km, 5km and 2km respectively; ESQRT equals 7.99 x 10 at this point. 
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MODIFIED SIMPSON AND ROBINSON MODEL 


Figure 11 
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Figure 12 
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C = 12.7 coulombs 


ELECTRIC FIELD MAP FOR 
MODIFIED MALAN MODEL 


Figure 13 
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VI , CONCLUSION 


In order to select the model which most, closely approximates a typica-1 
thundercloud, one must use more than the values of ESQRT which were 
generated. The modified models come closer to fitting the electric 

r 

field values measured at ground level by Wormell and the values Gish 
and Wait measured to 40,000 feet than the original models as shown in 
Figures 14 and 15, yet the final model should most nearly conform to 
all electric field measurements made in the vicinity of a thundercloud 
without lightning present. 

R. Markson [11] states that thunderclouds, not considering lightning 
flashes, can effect the electric field within 50km or more. In order 
to evaluate the ability o'f the modified models to meet this 
requirement, one can note the effect the various models have on the 
electric field at the farthest distance evaluated, i.e. 21km. Using, 
an expression for the approximate electric field given by Dolezalek [1], 
h being the height in km, 

E = 81.8 exp(-4.52h) + 38.6 exp(-.*375h) + 10.27 exp(-.121h) (9) 

the field values at 4km, 10km and 16km are 4.59 v/m, .640 v/m and 
.119 v/m respectively. The values at 21km horizontal distance from 
the charge centers at these heights are 16.53 v/m at 73° from the vertical, 
6.77 v/m at -12° from vertical and 4.24 v/m at -47° from vertical for the 
modified Simpson and Robinson model. The modified Gish and Wait model 
yi'elds 15.8 v/m at -0.1°, 5.18 v/m at -10° and 2.53 v/m at -36° and 



ELECTRIC FIELD STRENGTH (kv/m) 
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ELECTRIC FIELD STRENGTH AT GROUND LEVEL FOR THE 
MODIFIED MODELS 

Figure 15 
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34.42 'V/m at 20°, 45.2 v/m at -12° and 51.15 v/m at -56° for the 
modified Malan model. Since these values are for 21km and the field' 
is .to be effected within 50km or more, the largest resultant field 
values at the greatest angles must be chosen as being the closest 
fit. Here the modified Malan model is the best fin, followed by the 
modified Simpson and Robinson model, then the modified Gish and Wait 
model . 

D.R. Lane-Smith [12] states that "tropical thunderstorms have been 
observed to produce negative electric fields of 3000 v/m or higher." 

All three modified models qualify under this criterion, Simpson and 
Robinson having a negative field as high as 47,000 v/m at a distance 

I 

and height of 1km and 4 km, Gish and Wait producing 13,000 v/m at 1km 
and 4km height, and Malan producing 179,000 v/m at 1km and 10km height. 

C.G. Stergis, G.C. Kein and T. Kangas [13] in balloon flights over 
thunderclouds made the measurements shown in Figures 16 and 17 The 
electric field measurements at a horizontal distance of 1km and a 
height of 16km for the three modified models are 85.17 v/m for Simpson 
and Robinson, 33.05 v/m for Gish and Wait and 2886.61 v/m for Malan. 
Because the value of the field strength will decrease as the height 
increases to that of the readings in Figures 16 and 17, i.e. above 
21 km, the modified Malan model is the only one which could possibly 
satisfy these readings. 









37 


Height above the thundercloud 




500-1000 

1 ,000 

where E was measured, 






E n»x* v/m 




2,000 

2,500 

Number of measurements 

34 

29 

40 

24 

13 


DEPENDENCE OF E^ ON HEIGHT ABOVE THUNDERCLOUD 

Table 5 

(From Shvarts [14]) 


Similarly, Ya. M. Shvarts [14] gives the mean values of the maximum 
vertical component of the electric field at various distances above 
thunderclouds as shown in Table 5. In the modified Simpson and 
Robinson model the highest charge center is at 6km, the top of the 
cloud could therefore be considered to be somewhat below 8km. The 
vertical component of the electric field reading at 1km horizontal 
distance and the 8km height is 4821 v/m and, therefore, within the 
range of values in Table 5. The modified Gish and Wait model has its 
highest charge center also near 6km. In this case the vertical 
component of the electric field reading at 1km and 8km height is 
1822 v/m. This is somewhat lower than the values in Table 5. The 
uppermost charge center of the modified Malan model is at 10km. The 
top of the cloud could be considered to be near 13km. Therefore, using 
the reading at 1km and 14km height, the vertical component is 7239 v/m. 
This is slightly high yet still within the range shown by Shvarts. 
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When considering the electrical specifications for a typical thunderstorm 
as described in the last two sections, the modified Mai an model appears 
to be the best fit. The new electrical model of a cumulonimbus cloud 
has a charge of 20.1 coulombs at lukm, -36.8 coulombs at 5km and 
12.7 coulombs at 2km. 

The knowledge of the electric field' surrounding a typical thundercloud 
has some practical uses and some limitations. Low cost avionics such 
as the field measuring device as described by M.L. Hill [14] and later 
by R. Markson [11] would be able to detect a field mapping of a 
thundercloud and therefore warn the pilot to avoid such an area. Yet, 
this model is only of a typical thundercloud; any specific thundercloud 
could exhibit quite different electric field characteristics, particularly 
during lightning, which has not been accounted for in this model. A 
model which would include lightning would demand an exceedingly 
complex system of which this model is only a beginning. 
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APPENDIX A 

The following program, written for the HP2100 computer, is used to 
produce a field map for a system of three charge centers of values 
Q(1)»-Q{2) and Q(3) at the corresponding heights of V(l), Y(2) and 
Y(3). The output lists the X and Y locations, X field value, 

Y field value, total field and the angle in radians of the total 
field at the location (X,Y). X is incremented by 2000m from 1000m 
to 21000m;Y is incremented by 2000m from ground level to 16000m. 
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PROGRAM EMAP 
REAL NX,NY,KNY 
DIMENSION Y(6) , Q(6) 

20 READ (5,1) Y(l),Y(2),Y(3) f Q(l),Q(2),Q(3) 

1 FORMAT (6F10.2) 

IF ( Y ( 1 ) ) 22 , 22,21 

21 Y (4)=-Y (3) 

Y (5)— — Y (2) 

Y (6)=-Y (1 ) 

Q(4)=-Q(3) 

Q (5 )=-Q (2 ) 

Q(6)=-Q(l) 

READ (5,4) AT ,A2,A3‘,A4,A5,A6,A7 

4 FORMAT (7A4) 

WRITE (6,5) AT ,A2,A3,A4,A5 ,A6 ,A7 

5 FORMAT (1H1 ,A5X,7A4,//, 3X, "CHARGE CENTER 1 " ,5X , "CHARGE CENTER 2", 5 
IX , "CHARGE CENTER 3" ,/ ,2X , "HEIGHT" ,3X , "COULOMBS" , 3X , "HEIGHT" , 3X , "CO 
2UL0M8S" ,3X," HEIGHT" ,3X, "COULOMBS") 

WRITE (6,2) Y(1),Q(1),Y(2),Q(2),Y(3),Q(3) 

2 FORMAT (1H.6F10.2) 

WRITE (6,6) 

6 FORMAT (// ,3X,"X" ,5X,"Y" ,4X,"X FIELD", 3X,"Y FIELD" ,4X, "TOTAL ",5X," 
1 THETA") 

CONST = 4.*3.1416*8.854~(10.**(-12l) 

DO 10 NNY=0,1 6 ,2 
NY=NNY*1000. 

DO 10 NNX=1 0,2.1 0,20 
NX=NNX*100. 

EXSUM=0. 

EYSUM=0. 

DO 11 N=1 ,6,1 

EXSUM=EXSUM +Q(N)*MX/ (CONST* (NX* x 2+(NY-Y (N))- (NY-Y (N )))**! .5) 

1 1 EYSUM=EYSUM+Q (N )* (NY-Y (N ) )/ (CONST* (NX**2+(NY-Y (N)) x (NY-Y (N ) ) Hi . 5 
1) 

KNY=NY/1000. 

EYSUM=EYSUM-81 .8*EXP(-4.52*KNY)-38.*EXP(-.375*KNY)-1 0.27*EXP (- 121 
1*KNY) 

IF (EXSUM) 15,13,12 

12 THETA=ATAN(EYSUM/EXSUM) 

16 ETOTL=SQRT (EXSUM*EXSUM+EYSUM*EYSUM) 

14 WRITE (6,3) NX, NY, EXSUM, EYSUM,ETOTL, THETA 

3 FORMAT (216,3F10.2,F10.4) 

10 CONTINUE 

GO TO 20 

15 P 1=3 . 1 41 6 
ANGLE=SIGN(PI ,EYSUM) 

THETA=ANGLE+ATAN (EYSUM/EXSUM) 

GO TO 16 

!3 ANGLE=3. 1416/2. 

THETA=SIGN (ANGL E , EYSUM ) 

ETOTL=EYSUM 
GO TO 14 

22 END 





SIMPSON AND ROBINSON MODEL 


CHARGE CENTER 1 
HEIGHT COULOMBS 
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Y 
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4000 

1000 

6000 

3000 
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5000 
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7000 

6000 


CHARGE CENTER 2 
HEIGHT COULOMBS 

3000.00 -20.00 


X FIELD Y FIELD 


.00 

4066.83 

.00 

2559.72 

.00 

-881.20 

.00 

-1285.40 

.00 

-1054.19 

.00 

-797.06 

.00 

-605.92 

.00 

-474.12 

.00 

-383.90 

.00 

-321 .34 

.00 

-277.11 

-34552.17 

64956.74 

-7460.80- 

-123.11 

-787.46 

-1881.66 

245.10 

-1399.23 

306.39 

-967.06 

229.31 

-675.36 

157.24 

-483.80 

106.70 

-356.71 

73.30 

-270.68 

51 .36 

-211 .07 

36.75 

-168.81 

-^2327.79 

-97247.95 

-1222.37 

-13244.89 

1358.17 

-3666.89 

1077.23 

-1591.73 

714.07 

-902.53 

463.20 

-589.46 

303.79 

-415.05 

203.43 

-305.87 

139.51 

-232.83 

97.98 

-181.91 

70.38 

-145.35 

210423 44 

-15336.63 

17672.71 

-6032.24 

4772.50 

-2261 .52 

2042.89 

-1072.73 


CHARGE CENTER 3 
HEIGHT COULOMBS 

1500.00 4.00 


TOTAL 

THETA 

4066.83 

■ 1.5708 

2559.72 

1.5708 

-881 .20 

-1.5708 

-1285.40 

-1 5708 

-1054.19 

-1.5708 

-797.06 

-1 .5108 

-605.92 

-1.5708 

-474.12 

-1 .5708 

-383.90 

-1.5708 

-321 .34 

-1 .5708 

-277.11 

-1.5708 

73574.67 

2.0597 

7451 .82 

-3.1251 

2039.79 

-1 .9671 

1420.54 

-1.3974 

1014.43 

-1 .2640 

713.22 

-1.2435 

508.71 

-1.2566 

372.33 

-1 .2801 

280.43 

-1.3063 

217.23 

-1 .3321 

172.76 

-1 3564 

106060.39 

-1 .9813 

13301 .17 

-1 .6628 

3910.34 

-1 .2161 

1921.99 

-.9758 

1150.84 

-.9015 

749.68 

- . 9048 

514.35 

-.9390 

367.34 

-.9839 

271 .43 

-1 .0310 

206.62 

-1 .0767 

161.49 

-1.1199 

210981.59 

-.0728 

18673.86 

-.3289 

5281.21 

- 4425 

2307.41 

-.4835 
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X 

Y 

X FIELD 

Y FI-ELD 

TOTAL 

THETA 

9000 

6000 

1093.99 

-631 .81 

1263.32 

-.5237 

11000 

6000 

655.36 

-427.21 

782.31 

-.5777 

13000 

6000 

418.42 

-311.83 

521 .84 

-.6405 

15000 

6000 

278.75 

-237.61 

366.28 

-.7059 

17000 

6000 

191.85 

-186.06 

267.25 

-.7701 

19000 

6000 

135.65 

-148.64 

201 .24 

-.8310 

21000 

6000 

98.18 

-120.74 

1 55-. 62 

-.8881 

1000 

8000 

18077.88 

32606.24 

37282.39 

1.0646 

3000 

8000 

11414.12 

5261 .37 

12568.38 

.4319 ~ 

5000 

8000 

4725.81 • 

563.86 

4759.33 

.1188' 

7000 

8000 

2246.16 

-134.77 

2250.20 

-.0599 

9000 

8000 

1229.38 

-226.77 

1250.12 

-.1824 

11000 

8000 

.742.58 

-215.60 

773.25 

-.2826 

13000 

8000 

478.86 

-187.20 

514.15 

- . 3' 7 27 

15000 

8000 

323.02 

-158.48 

359.85 

-.^564 

17000 

8000 

225.34 

-133.35 

261 .84 

-.5343 

19000 

8000 

161.45 

-112.08 

196.54 

-.6068 

21000 

8000 

118.28 

-94.50 

151.40 

-.6741 

1000 

10000 

2631 .71 

9182.00 

9551 .70 

1.2917 

3000 

10000 

4116.18 

4404.92 

6028.78 

.8193 

5000 

10000 

2896.25 

1547.05 

3283.54 

.4906 

7000 

10000 

1777.12 

483.03 

1841.60 

.2654 

9000 

10000 

1102.76 

111.90 

1108.42 

.1011 

1100.0 

10000 

712.94 

-20.71 

713.24 

-.0290 

13000 

10000 

479.75 

-67.37 

484.45 

-.1395 

15000 

10000 

333.74 

-80.73 

343.37 

-.2373 

17000 

10000 

238.55 

-80.62' 

251 .80 

-.3259 

19000 

10000 

174.38 

-75.18 

189.90 

-.4071 

21000 

10000 

129.96 

-67.94 

' 146.64 

-.4817 

1000 

12000 

748.58 

3825.96 

3898.50 

1.3776 

3000 

12000 

1598.30 

2610.45 

3060.89 

1.0214 

5000 

12000 

1556.23 

1395.26 

2090.12 

.7309 

7000 

12000 

1199.36 

662.58 

1370.21 

.5047 

9000 

12000 

860.01 

291.66 

908.12 

.3270 

11000 

12000 

609.46 

112.35 

619.73 

.1823 

13000 

12000 

435.98 

26.07 

436.76 

.0597 

15000 

12000 

316.79 

-15.06 

317.14 

-.0475 

17000 

12000 

234.00 

-33.81 

236.43 

-.1435 

19000 

12000 

175.57 

-41.23 

180 35 

-.2307 

21000 

12000 

133.65 

-42.89 

140.36 

-.3106 

1000 

14000 

296.39 

1984.18 

2006.20 

1.4225 

3000 

14000 

723.55 

1573.04 

1731 .51 

1 .1396 

5000 

14000 

845.55 

1044.32 

1343.71 

.8902 

7000 

14000 

767.55 

620.30 

986.87 

.6797 

9000 

14000 

623.99 

344.12 

712.58 

.5040 

11000 

14000 

484.45 

179.98 

516.81 

.3557 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

13000 

14000 

370.29 

85.95 

380.14 

.2281 

15000 

14000 

282.64 

32.89 

284.55 

.1159 

17000 

14000 

216.82 

3.36 

216.84 

.0155 

19000 

• 14000 

167. ,62 

-12.68 

168.10 

-.0755 

21000 

14000 

130.74 

-20.94 

132.4T 

-.1588 

1000 

16000 

142.35 

1173.79 

1182.39 

1.4501 

3000 

16000 

372.18 

1004.45 

1071.18 

1 .2160 

5000 

16000 

483.71 

754.30 

896.07 

1 .0006 

7000 

16000 

490.83 

515.07 

711 .48 

' .8095 

9000 

16000 

440.14 

329.57 

549.85 

.6427 - 

11000 

16000 

370.05 

200.88 

421 .06 

.4973 

13000 

16000 

301.17 

116.66 

322.98 

.3696 

15000 

16000 

241.51 

63.21 

2*9.64 

.2560 

17000 

16000 

192.67 

29.86 

194.97 

.1537 

19000 

1.6000 

153.74 

9.36 

154.03 

.0608 

21000 

16000 

123.08 

-3.00 

123.12 

-.0244 
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GISH AND WAIT MODEL 


CHARGE 

CENTER 1 

CHARGE 

CENTER 2 

CHARGE 

HEIGHT 

COULOMBS 

HEIGHT 

COULOMBS 

HEIGHT 

6096.00 

39.00 

3048.00 

-39.00 

.00 


X 

Y 

X FIELO 

Y FIELD 

. TOTAL 

1000 

0 

-.00 

46473.19 

46473.19 

3000 

0 

.00 

13560.42 

13560.42 

5000 

0 

-.00 

1791 .56 

1791 .56 

7000 

0 

.00 

-672.18 

-672.18 

9000 

0 

.00 

-966.69 

-966.69 

11000 

0 

.00 

-841 .81 

-841 .81 

13000 

0 

.00 

-676.21 

-676.21 

15000 

0 

.00 

-541.02 

-541.02 

17000 

0 

.00 

-440.81 

-440.81 

19000 

0 

-.00 

-368.09 

368.09 

21000 

0 

.00 

-315.17 

-315.17 

1000 

2000 

-108718.89 

109432.34 

154257.06 

3000 

2000 

-21174.68 

4780.39 

21707.59 

5000 

2000 

-3802.59 

-948 . 36 

3919.07 

7000 

2000 

-504.58 

-1243.51 

1341 .98 

9000 

2000 

114.77 

-1006.48 

1013.00 

11000 

2000 

187.85 

-757.27 

780.22 

13000 

2000 

157.06 

-565.07 

586.50 

15000 

2000 

116.66 

-426.54 

442.21 

17000 

2000 

84.33 

-328.00 

338.67 

19000 

2000 

61 .00 

-257.49 

264.62 

21000 

2000 

44.61 

-206.39 

211.15 

1000 

4000 

-104554.09 

-181997.62 

209892.09 

3000 

4000 

-10832.73 

-23238.61 

25639.44 

5000 

4000 

-803.15 

-5802.10 

5857.43 

7000 

4000 

509.48 

-2232.55 

2289.94 

9000 

4000 

578.85 

-1172.51 

1307.61 

11000 

4000 

452.19 

-741.18 

868 24 

13000 

4000 

325.58 

-516 84 

610.84 

15000 

4000 

230.37 

-380.52 

444.82 

17000 

4000 

163.69 

-290.07 

333.07 

19000 

4000 

117.79 

-226.91 

255.66 

21000 

4000 

86.10 

-181 .33 

200.73 

1000 

6000 

334420.31 

-65541.23 

340782.31 

3000 

6000 

25453.36 

-13660.82 

28887.57 

5000 

6000 

5864.96 

-4583.87 

7443.76 


CENTER 3 
COULOMBS 
.00 


THETA - - 


1 .5708 
1.5708 
1.5708 
-1 .5708 
-1.5708 
-1 .5708 
-1.5708 
-1 .5078 
-1 .5708 
-1.5708 
-1.5708 
2.3529 
2.9196 
-2.8972 
-1 .9563 
-1 4573 
-1.3276 
-1 .2997 
-1.3038 
-1 .3191 
-1 .3382 
-1 .3579 
-2.0922 
-2.0070 
-1 .7084 
-1 .3464 
-1 .1122 
-1 .0230 
-1 .0087 
-1.0264 
-1 .0570 
-1 .0920 
-1.1275 
-.1935 
-.4926 
-.6634 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

7000 

6000 

2295.37 

-1901 .85 

2980.90 

-.6919 

9000 

6000 

1211.50 

-984.07 

1560.81 

-.6822 

11000 

6000 

738.83 

-607.15 

956.30 

-.6879 

13000 

6000 

483.46 

-420.12 

640.50 

-.7154 

15000 

6000 

329.47 

-311.14 

453.16 

-.7568 

17000 

6000 

231 .05 

-240.02 

333.16 

-.8044 

19000 

6000 

165.83 

-190.20 

252.34 

-.8537 

21000 

6000 

121 .47 

-153.73 

195.93 

-.9021 

1000 

8000 

32652.30 

54711 .66 

63714.52 

1.0327 

3000 

8000 

18372.59 

6858.57 

19611.02 

.3573 - 

5000 

8000 

6873.42 

65.88 

6873.7& 

.0096 

7000 

8000 

3002.36 

-547.50 

3051 .88 

-.1804 

9000 

8000 

1557 87 

-467.76 

1626.58 

-.2917 

11000 

8000 

917.97 

-357.29 

985.05 

-.3712 

13000 

8000 

587.85 

-277.19 

649.92 

-.4406 

15000 

8000 

397.35 

-220.69 

454.52 

-.5070 

17000 

8000 

278.74 

-179.23 

331 .39 

-.5714 

19000 

8000 

201 .04 

-147.64 

249.43 

-.6334 

21000 

8000 

418.25 

-122.96 

192.61 

-.6924 

1000 

10000 

4416.47 

14565:66 

15220.50 

1.2764 

3000 

10000 

6586.64 

6468.16 

9231 .52 

.7763 

5000 

10000 

4350.91 

1974.11 

4777 82 

.4259 

7000 

10000 

2512.04 

481.38 

2557.75 

.1893 

9000 

10000 

1485.50 

37.40 

1485.97 

.0252 

11000 

10000 

930.13 

-90.61 

934.53 

-.0971 

13000 

10000 

614.63 

-122.40 

626.70 

-.1966 

15000 

10000 

423.93 

-123.41 

441 .53 

-.2833 

17000 

10000 

302.16 

-114.40 

323.09 

-.3619 

19000 

10000 

220.97 

-102.57 

243 62 

-.4346 

21000 

10000 

165.01 

-90.60 

188.25 

-.5021 

1000 

12000 

1196.20 

5814.28 

5936 05 

1.3679 

3000 

12000 

2491.73 

3826.18 

4566.00 

.9936 

5000 

12000 

2333.49 

1919.59 

3021 .59 

.6884 

7000 

12000 

1724.09 

839.38 

1917 56 

.4531 

9000 

12000 

1190.52 

332.74 

1236.15 

.2725 

11000 

12000 

819.34 

107.33 

826.34 

.1303 

13000 

12000 

574.33 

7.58 

574.38 

.0132 

15000 

12000 

412.00 

-36.05 

413.58 

-.0873 

17000 

12000 

302.12 

-53.86 

306.89 

-.1764 

19000 

12000 

225.87 

-59.41 

233.56 

-.2572 

21000 

12000 

171 .73 

-59.08 

181 .61 

-.3313 

1000 

Uooo 

456.98 

2925.56 

2961 .03 

1 .4158 

3000 

14000 

1099.50 

2273.60 

2525.50 

1 .1204 

5000 

14000 

1253.13 

1456.99 

1921 .76 

.8605 

7000 

14000 

1104.81 

827.48 

1380 34 

.6429 

9000 

14000 

873.05 

436.67 

976.17 

.4638 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

11000 

14000 

661 .57 

215.99 

695.94 

.3156 

13000 

14000 

496.24 

95.67 

505.37 

.1905' 

15000 

14000 

373.66 

30.78 

374.92 

.0822 

17000 

14000 

284 ;04 

-3.87 

284.06 

-.0136 

19000 

14000 

218.33 

-21.90 

219.43 

-.1000 

21000 

14000 

169.73 

-30.69 

172.48 

-.1789 

1000 

16000 

213.59 

1693.36 

1706.78 

1 .4453 

3000 

16000 

553.23 

1431.78 

1534.94 

1.2021 

5000 

16000 

707.20 

1052.13 

1267.72 

.9790 

7000 

16000 

703.10 

698.73 

991.25 

.7823 

9000 

16000 

617.32 

433.55 

754.35 

.6123 

11000 

16000 

509 07 

256.01 

569.82 

.4660 

13000 

16000- 

407.65 

143.80 

432.27 

.3391 

15000 

16000 

322.76 

74.81 

331 .31 

.2278 

17000 

16000 

255.06 

32.95 

257.18 

.1285 

19000 

16000 

202.18 

7.85 

202.34 

.0388 

21000 

16000 

161.14 

-6.94 

161 .28 

-.0431 
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MAL AN MODEL 


CHARGE CENTER 1 

CHARGE 

CENTER 2 

CHARGE 

HEIGHT 

COULOMBS 

HEIGHT 

COULOMBS 

HEIGHT 

10000.00 

40.00 

5000.00 

-40.00 

2000.00 


X 

Y 

X FIELD 

Y FIELD 

TOTAL 

1000 

0 

.00 

-12251 .73 

12251 .73 

3000 

0 

-.00 

4015.52 

4015.52 

5000 

0 

.00 

2591 .42 

2591.52 

7000 

0 

-.00 

632.48 

632.48 

9000 

0 

.00 

-247.33 

-247.33 

11000 

0 

.00 

-537.95 

-537.95 

13000 

0 

.00 

-587.30 

-587.30 

15000 

0 

.00 

-551.52 

-551 .52 

17000 

0 

.00 

-492.79 

-492.79 

19000 

0 

.00 

-433.99 

-433.99 

21000 

0 

.00 

-382.39 

-382.39 

1000 

2000 

78723.22 

28109.79 

83591.31 

3000 

2000 

-2692.49 

10027.17 

10382.37 

5000 

2000 

-3037.11 

2608.98 

4003.85 

7000 

2000 

-1313.15 

320.30 

1351 .65 

9000 

2000 

-449.35 

-366.78 

580.04 

11000 

2000 

-102.73 

-530.26 

540.12 

13000 

2000 

21.59 

-520.62 

521 .07 

15000 

2000 

58.41 

-459.97 

463.67 

17000 

2000 

62.82 

-391.33 

396.34 

19000 

2000 

56.43 

-328.87 

333.68 

21000 

2000 

47.45 

-276.16 

280.20 

1000 

4000 

-117514.45 

133725.22 

178022.72 

3000 

4000 

-24778.07 

8331 .59 

26141 .32 

5000 

4000 

-6751 .43 

-378.73 

6762.04 

7000 

4000 

-2037.31 

-1102.32 

2316.41 

9000 

4000 

-540.37 

-954.34 

1096.70 

11000 

4000 

-42.81 

-754.78 

756.00 

13000 

4000 

111.29 

-595.46 

605.77 

15000 

4000 

144.15 

-474.77 

496.17 

17000 

4000 

135.96 

-382.97 

406.39 

19000 

4000 

116.25 

-312.38 

333.31 

21000 

4000 

95.57 

-257.55 

274.71 

1000 

6000 

-120686.14 

-142335.00 

186612.97 

3000 

6000 

-23275.31 . 

-19824.02 

30573.39 

5000 

6000 

-4897.77 

-6664.76 

8270.87 


CENTER 3 
COULOMBS 
10.00 


THETA 


-1.5708 

1 . 5708 

1 . 5708 

1.5708 
-1.5708 
-1.5708 
-1.5708 
-1.5708 
-1.5708 
-1.5708 
-1.5708 
.3430 
1.8331 
2.4319 
2.9024 

-2.4570 
-1.7622 
-1.5294 
-1 .4445 
-1.4116 
-1.4009 
-1.4006 
2.2918 
2.8172 
-3.0856 
-2.6456 
-2.0860 
-1.6275 
-1.3860 
-1.2760 
-1.2297 
-1.2145 
-1.2155 
-2.2741 
-2.4361 
-2.2046 
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X 

Y 

X FIELD 

7000 

6000 

-976.07 

9000 

6000 

16.83 

11000 

6000 

265.29 

13000 

6000 

297.59 

15000 

6000 

266.77 

17000 

6000 

222.12 

19000 

6000 

179.59 

21000 

6000- 

143.60 

1000 

8000 

21198.36 

3000 

8000 

9820.07 

5000 

8000 

3454.41 

7000 

8000 

1714.37 

9000 

8000 

1067.79 

11000 

8000 

742.80 

13000 

8000 

544.33 

15000 

8000 

409.54 

17000 

8000 

313.02 

19000 

8000 

241 .97 

21000 

8000 

188.84 

1000 

10000 

356977.94 

3000 

10000 

34966.09 

5000 

10000 

9876.40 

7000 

10000 

3963.03 

9000 

10000 

1994.93 

11000 

10000 

1175.08 

13000 

10000 

768.44 

15000 

10000 

537.57 

17000 

10000 

392.89 ' 

19000 

10000 

295.74 

21000 

10000 

227.34 

1000 

12000 

31233.70 

3000 

12000 

20824.47 

5000 

12000 

9037.53 

7000 

12000 

4310.54 

9000 

12000 

2316.31 

11000 

12000 

1381.88 

13000 

12000 

895.93 

15000 

12000 

618.29 

17000 

12000 

446.61 

19000 

12000 

333.48 

21000 

12000 

255.22 

1000 • 

14000 

4700.79 

3000 

14000 

7520.32 


Y FIELD 

TOTAL 

THETA 

-2978.22 

3134.08 

-1.8875 

-1586.22 

1586.31 

-1.5602 

-969.45 

1005.09 

-1.3037 

-657.70 

721.90 

-1.1459 

-480.23 

549.35 

-1.0637 

-368.41 

430.91 

-1.0282 

-292.18 

342.97 

-1.0197 

-237.19 

277.27 

-1.0264 

-95907.80 

98222.58 

-1.3533 

-27568.87 

29265.62 

-1.2286 

-8824.13 

9475.19 

-1.1977 

-3563.27 

3954.23 

-1.1224 

-1741.09 

2042.45 

-1.0207 

-990.34 

1237.95 

-.9273 

-634.98 

836.36 

-.8621 

-445.56 

604.19 

-.8275 

-333.47 

457.35 

-.8170 

-260.99 

355.90 

-.8232 

-210.66 

282.91 

-.8399 

-12116.47 

357183.50 

-.0339 

-7850.73 

35836.60 

-.2209 

-4179.14 

10724.20 

- . 4003 

-2198.17 .- 

4531 .83 

-.5064 

.-1232.97 

2345.20 

-.5536 

-753.90 

1396.13 

-.5704 

-502.01 

917.89 

-.5787 

-359.77 

646.85 

- . 5898 

-273.12 

478.49 

-.4076 

-216.30 

366.^0 

-.6315 

-176.57 

287.85 

-.6604 

58117 42 

65978.62 

1.0777 

10466.50 

23306.79 

.4657 

1317.70 

9133.19 

.1448 

-266.50 

4316.48 

-.0525 

-417.18 

2353.58 ' 

-.1782 

-368.39 

1430.15 

-.2605 

-295.87 

943.52 

-.3190 

-236.96 

662.14 

-.3650 

-193.30 

486.64 

-.4085 

-160.96 

370.29 

-.4497 

-136.36 

289.36 

-.4907 

16794.50 

17439.98 

1.2979 

8299.14 

11199.60 

.8346 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

5000 

14000 

5424.10 

3011.47 

6204.01 

.5068 

7000 

14000 

3362.94 

961 .95 

3497.82 

.2786 

9000 

14000 

2082.60 

243.79 

2096.82 

.1165 

11000 

14000 

1341 .62 

-2.14 

1341 .62 

-.0016 

13000 

14000 

905.85 

-83.00 

909.65 

-.0914 

15000 

14000 

638.60 

-105.42 

647.24 

-.1636 

17000 

14000 

466.70 

-107.03 

478.82 

-.2254 

19000 

14000 

351.07 

-101.35 

365.41 

-.2810 

21000 

14000 

270.20 

-93.40 

285.89 

-.3328 ' 

1000 

16000 

1366.02 

7108.68 

7 238.74 

1.3809 

3000 

16000 

2945.49 

4908.01 

5724.03 

1.0303 

5000 

16000 

2903.74 

2666.07 

3942.04 

.7428 

7000 

16000 

2255.76 

1287.59 

2597.37 

.5187 

9000 

16000 

1621 .80 

584 09 

1723.78 

.3457 

11000 

16000 

1150.72 

246.31 

1176.79 

.2109 

13000 

16000 

826.67 

85.77 

831.11 

.1034 

15000 

16000 

606.64 

9.07 

606.71 

.0150 

17000 

16000 

455.31 

- 27.59 

456.15 

-.0605 

19000 

16000 

348.89 

-44.61 

351 .73 

-.1272 

21000 

16000 

272.20 

-51.70 

277.07 

-.1877 
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MODIFIED SIMPSON AND ROBINSON MODEL 


CHARGE CENTER 1 
HEIGHT COULOMBS 

6000.00 4.60 


X 

Y 

1000 

0 

3000 

0 

5000 

0 

7000 

0 

9000 

0 

11000 

0 

13000 

0 

15000 

0 

17000 

0 

19000 

0 

21000 

0 

1000 

2000 

3000 

2000 

5000 

2000 

7000 

2000 

9000 

2000 

11000 

2000 

13000 

2000 

15000 

2000 

17000 

2000 

19000 

2000 

21000 

2000 

1000 

4000 

3000 ■ 

4000 

5000 

4000 

7000 

4000 

9000 

4000 

11000 

4000 

13000 

4000 

15000 

4000 

17000 

4000 

19000 

4000 

21000 

4000 

1000 

6000 

3000 

6000 

5000 

6000 

7000 

6000 

9000 

6000 


CHARGE CENTER 2 
HEIGHT COULOMBS 

3000.00 -11.30 


X FIELD Y FIELD 


.00 

-10538.48 

.00 

1939.86 

-.00 

770.68 

-.00 

177.66 

.00 

-20.31 

-.00 

-88.79 

-.00 

-114.33 

-.00 

-124.48 

.00 

-218.68 

.00 

-130.42 

-.00 

-131.11 

2650.72 

52024.13 

-3226.18 

2930.23 

-918.40 

537.89 

-279.52 

153.89 

-93.95 

41 .60 

-34.56 

.63 

-13.63 . 

-15.77 

-5.63 

-22.64 

-2.35 

-25.57 

-.96 

-26.79 

-.34 

-27.25 

■29521.43 

-36532.12 

-4359.65 

-2692.04 

974.39 

-440.09 

-304.88 

-102.55 

-111.29 

-34.15 

-44.17 

-19.65 

-18.43 

-16.91 

-7.89 

-16.60 

-3.35 

-16.64 

-1.32 

-16.62 

-.41 

-16.52 

38668.84 

-7161.18 

1607.50 

-2427.45 

-21.39 

-690.32 

-90.12 

-220.81 

-53.42 

-82.80 


CHARGE CENTER 3 
HEIGHT COULOMBS 

1500.00 5.97 


TOTAL 

THETA 

-10538.48 

-1.5708 

1939.85 

1.5708 

770.68 

1 .5708 

177.66 

1 .5708 

-20.31 

-1 .5708 

88.79 

-1.5708 

114.33 

-1.5708 

• 124.48 

-1.5708 

-128.68 

-1.5708 

-130.42 

-1 .5708 

131.11 

-1.5708 

52091 .62 

1.5199 

4358.27 

2.4042 

1064.33 

2.6118 

319.08 

2.6383 

102.75 

2.7247 

34.56 

3.1235 

20.84 

-2.2837 

23.33 

-1 .8145 

25.68 

-1.6630 

26.81 

-1.6066 

27.25 

-1.5832 

46969.25 

-2.2505 

5123.83 

-2.5884 

1069.16 

-2.7174 

321.66 

-2.8171 

116.41 

-2.8439 

48.34 

-2.7231 

25.02 

-2.3993 

18.38 

-2.0143 

16.97 

-1 .7693 

16.67 

-1.6502 

16.53 

-1.5956 

39326.34 

-.1831 

2911.46 

-.9859 

690.65 

-1 .6018 

238.49 

-1.9583 

98.54' 

-2.1438 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

11000 

6000 

-26.85 

-37.86 

46.42 

-2.1877 

- 13000 

6000 

-12.78 

-21.92 

25.37 

-2.0985 

15000 

6000 

-5.82 

-15.78 

16.81 

-1.9240 

17000 

6000 

-2.44 

-13.17 

13.39 

-1.7544 

19000 

6000 

-.83 

-11.92 

11.95 

-1.6404 

21000 

6000 

-.08 

-11.23 

11.23 

-1.5778 

1000 

8000 

3119.30 

4820.74 

5741.91 

.9965 

3000 

8000 

1548.54 

188.47 

1559.96 

.1211 . 

5000 

8000 

381.93 

-235.11 

448.49 

-.5518 

7000 

8000 

80.84 

-145.22 

167.08 

-1.0658' 

9000 

8000 

12.37 

-74.57 

75.59 

-1.4065 

11000 

8000 

-.95 

-39.09 

39.10 

-1.5951 

13000 

8000 

-2.06 

-22.69 

22.79 

-1.6615 

15000 

8000 

-1.10 

-15.06 

15.10 

-1.6439 

17000 

8000 

-.20 

-11.36 

11.36 

-1.5885 

19000 

8000 

.34 

-9.46 

9.47 

-1.5350 

21000 

8000 

.60 

-8.41 

8.^3 

-1.5000 

1000 

10000 

389.07 

1106.02 

1172.46 

1.2325 

3000 

10000 

525.92 

367.35 

641.52 

.6097 

5000 

10000 

277.42 

19.99 

278.14 

.0719 

7000 

10000 

113.39 

-46.45 

122 53 

-.3888 

9000 

10000 

43.38 

-41.14 

59.79 

-.7589 

11000 

10000 

16 93 

-27.92 

32.65 

-1.0257 

13000 

10000 

7.27 

-18.44 

19.82 

-1.1950 

15000 

10000 

3.73 

-12.78 

13.31 

-1 .2868 

17000 

10000 

2.36 

-9.55 

9.84 

-1.3285 

19000 

10000 

1.76 

-7.71 

7.90 

-1.3468 

21000 

10000 

1 .43 

-6.62 

6.77 

-1.3583 

1000 

12000 

93.86 

377.89 

389.37 

1.3274 

3000 

12000 

181.48 

210.26 

277.75 

.8587 

5000 

12000 

145.99 

66.95 

160.61 

.4300 

7000 

12000 

86.33 

3.99 

86.42 

.0462 

9000 

12000 

45.34 

-13.40 

47.28 

-.2873 

11000 

12000 

23.26 

-14.64 

27.48 

-.5617 

13000 

12000 

12.34 

-12.05 ' 

17.24 

-.7734 

15000 

12000 

7.04 

-9.38 

11.72 

-.9270 

17000 

12000 

4.40 

-7.41 

8.61 

-1.0351 

19000 

12000 

3.00 

-6.08 

6.78 

-1.1119 

21000 

12000 

2.21 

-5.22 

5.66 

-1.1707 

1000 

14000 

31.83 

164.82 

167.86 

1 . 3800 

3000 

14000 

72.40 

114.53 

135.50 

1.0071 

5000 

14000 

73.86 

56.83 

93.20 

.6558 

7000 

14000 

55.56 

19.21 

58.78 

.3328 

9000 

14000 

36.16 

1.52 

36.20 

.0419 

11000 

14000 

22.20 

-4.80 

22.71 

-.2130 
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, X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

13000 

14000 

13.54 

-6.19 

14.89 

-.4291 

15000 

14000 

8.48 

-5.88 

10.32 

-.6067 

17000 

14000 

5.56 

-5.18 

7.60 

- .7^99 

19000 

14000 

3.84 

-4.51 

5.92 

-.8656 

21000 

14000 

2.79 

-3.98 

4.86 

-.9603 

1000 

16000 

13.30 

84.13 

85.17 

1 .4140 

3000 

16000 

33.00 

65.71 

73.53 

1.1054 

5000 

16000 

38.80 

40.76 

56.27 

.8100 

7000 

16000 

34.25 

20.22 

39.77 

.5332 

9000 

16000 

26.01 

7.43 

27.05 

.2784 

11000 

16000 

18.29 

.87 

18.31 

.0476 

13000 

16000 

12.47 

-1.98 

12.63 

-.1578 

15000 

16000 

8.50 

-2.99 

9.01 

-.3378 

17000 

16000 

5.90 

-3.18 

6.71 

-.4942 

19000 

16000 

4.22 

-3.07 

5.22 

-.6293 

21000 

16000 

3.11 

-2.88 

4.24 

-.7463 
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MODIFIED GTSH AND WAIT MODEL 


CHARGE CENTER 1 
HEIGHT COULOMBS 

6096.00 1.52 


CHARGE CFNTER 2 
HEIGHT COULOMBS 
3048.00 -2.64 


CHARGE CENTER 3 
HEIGHT COULOMBS 

.00 .00 


X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA. . 

1000 

0 

.00 

3545.22 

3545.22 

1.57C8 

3000 

0 

-.00 

1187.98 

1187.98 

1.5708 

5000 

0 

.00 

250.43 

250.43 

1.5708 

7000 

0 

.00 

-13.31 

-13.31 

-1.5708 

9000 

0 

-.00 

-91 .15 

91.15 

-1 .5708 

11000 

0 

.00 

-116.55 

-116.55 

-1.5708 

13000 

0 

.00 

-125.58 

' -125.58 

-1.5708 

15000 

0 

.00 

-128.98 

-128.98 

-1.5708 

17000 

0 

.00 

-130.27 

-130.27 

-1 .5708 

19000 

0 

.00 

-130.73 

-130 73 

-1.5708 

21000 

0 

.00 

-130.86 

-130.86 

-1.5708 

1000 

2000 

-7475.19 

8083.68 

11010 19 

2.3171 

3000 

2000 

-1617.19 

741.00 

1778.87 

2.7119 

5000 

2000 

-385.38 

158.82 

416.83 

2.7507 

7000 

2000 

-108.76 

35.33 

114.35 

2.8275 

9000 

2000 

-34.87 

-3.81 

35.07 

-3.0328 

11000 

2000 

-12.29 

-18.04 

21.83 

-2.1690 

13000 

2000 

-4.62 

-23.54 

23.99 

-1 .7645 

15000 

2000 

-1.78 

-25.71 

25.77 

-1 .6399 

17000 

2000 

-.66 

-26.54 

26.55 

-1.5958 

19000 

2000 

-.21 

-26.82 

26.83 

-1.5785 

21000 

2000 

-.02 

-26.87 

26.87 

-1 .5716 

1000 

4000 

-7871.57 

-10551.76 

13164.40 

-2.2117 

3000 

4000 

-1323.56 

-1069.44 

1701 .62 

-2.4620 

5000 

4000 

-335.01 

-203.10 

391.77 

-2.5966 

7000 

4000 

-108.12 

-56.04 

121.78 

-2.6634 

9000 

4000 

-39.01 

-25.36 

46.53 

-2.5651 

11000 ' 

4000 

-14.98 

-18.42 

23.74 

-2.2536 

13000 

4000 

-5.93 

-16.76 

17.78 

-1.9109 

15000 

4000 

-2.32 

-16.30 

16.46 

-1.7123 

17000 

4000 

-.83 

-16.09 

16.12 

-1.6224 

19000 

4000 

-.21 

-15.94 

15.94 

-1 . 5838 

21000 

4000 

.05 

-15.80 

15.80 

-1.5677 

1000 

6000 

12714.70 

-3423.76 

13167.60 

-.2630 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL ' 

THETA 

3000 

6000 

621 .85 

-834.46 

1040.68 

-.9304 

5000 

6000 

17.03 

-256.63 

257.19 

-1.5045 

7000 

6000 

-24.18 

-89.68 

92.89 

-1.8342 

9000 

6000 

-15.81 

-38.13 

41 .28 

-1.9639 

nooj 

6000 

-7.83 

-20.88 

22.30 

-1.9297 

13000 

6000 

-3.46 

-14.62 

15.02 

-T.8032 

15000 

6000 

-1.33 

-12.12 

12.20 

-1.6801 

17000 

6000 

-.35 

-11.01 

11.02 

-1.6024 

19000 

6000 

.08 

-10.44 

10.45 

-1.5631 

21000 

6000 

.25 

-10.11 

10.12 

-1.5*62 • 

1000 

8000 

1201 .91 

1821 .65 

2182.42 

.9876 

3000 

8000 

580.59 

79.03 

585.95 

.1353 

5000 

8000 

151.69 

-83.67 

173.24 

-.5041 

7000 

8000 

36.74 

-56.27 

67.20 

-.9924 

9000 

8000 

8.45 

-31 .35 

32.47 

-1.3076 

11000 

8000 

1.92 

-18.54 

18.64 

-1.4675 

13000 

8000 

.66 

-12.49 

12.51 

-1.5183 

15000 

8000 

.53 

-9.61 

9.63 

-1.5156 

17000 

8000 

.57 

-8.19 

8.21 

-1.5011 

19000 

8000 

.58 

-7.43 

7.45 

-1 .4925 

21000 

8000 

.55 

-7.00 

7.02 

-1.4919 

1000 

10000 

147.57 

420.51 

445.66 

1.2333 

3000 

10000 

199.73 

141 .81 

244.95 

.6174 

5000 

10000 

107.87 

9.84 

108.32 

.0909 

7000 

10000 

46.23 

-17.46 

49.42 

-.3612 

9000 

10000 

19.05 

-16.94 

25.49 

-.7270 

11000 

10000 

8.28 

-12.67 

1 5.. 1 3 

-.9921 

13000 

10000 

4.05 

-9.39 

10.23 

-1.1635 

15000 

10000 

2.32 

-7.38 

7.74 

-1.2658 

17000 

10000 

1.54 

-6.21 

6.40 

-1.3272 

19000 

10000 

1.14 

-5.52 

5.64 

-1.3677 

21000 

10000 

.89 

-5.11 

5.18 

-1.3982 

1000 

12000 

35.78 

144.76 

149.11 

1.3285 

3000 

12000 

69.59 

81 .26 

106.99 

.8626 

5000 

12000 

56.88 

26.50 

62.75 

.4360 

7000 

12000 

34.59 

1.73 

34.64 

.0499 

9000 

12000 

18.92 

-5.70 

19.76 

-.2928 

11000 

12000 

10.21 

-6.73 

12.23 

-.5828 

13000 

12000 

5.75 

-6.07 

8.36 

-.8129 

15000 

12000 

3.47 

-5.24 

6.28 

-.9853 

17000 

12000 

2.27 

-4.58 

5.11 

-1.1104 

19000 

12000 

1.59 

-4.12 

4.42 

-1.2015 

21000 

12000 

1.18 

-3.81 

3.99 

-1.2695 

1000 

14000 

12.26 

63.59 

64.76 

1 . 3803 

3000 

V4000 

28.05 

44.40 

52.52 

1.0074 

5000 

14000 

28.97 

22.21 

36.51 

.6540 
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X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

7000 

14000 

22.23 

7.46 

23.45 

.3237 

9000 

14000 

14.87 

.26 

14.87 

.0178 

11000 

14000 

9.44 

-2.52 

9.77 

-.2606 

13000 

14000 

5.98 

-3.31 

6.83 

-.5054 

15000 

14000 

3.89 

3.36 

5.14 

-.7120 

17000 

14000 

2.64 

-3.19 

4.14 

-.8794 

19000 

14000 

1.87 

-2.99 

3.53 

-1.0117 

21000 

14000 

,1 .39 

-2.83 

3.15 

-1.1152 

1000 

16000 

5.18 

32.64 

33.05 

1.4133 

3000 

16000 

12.92 

25.55 

28.63 

1.1025 

5000 

16000 

15.35 

15.88 

22.08 

.8024' 

7000 

16000 

13.75 

7.81 

15.81 

.5162 

9000 

16000 

10.66 

2.66 

10.98 

.2445 

11000 

16000 

7.68 

-.09 

7.68 

-.0118 

13000 

16000 

5.38 

-1.37 

5.55 

-.2500 

15000 

16000 

3.78 

-1.90 

4.22 

-.4653 

17000 

16000 

2.69 

-2.06 

3.39 

-.6541 

19000 

16000 

1.97 

-2.09 

2.87 

-.8137 

21000 

16000 

1.48 

-2.05 

2.53 

-.9455 
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MODIFIED MALAN MODEL 


.CHARGE CENTER 1 CHARGE 
HEIGHT COULOMBS HEIGHT 
10000.00 20.10 5000.00 


X 

Y 

X FIELD 

1000 

0 

.00 

3000 

0 

-.00 

5000 

0 

-.00 

7000 

0 

-.00 

9000 

0 

.00 

11000 

0 

-.00 

13000 

0 

-.00 

15000 

0 

-.00 

17000 

0 

.00 

19000 

0 

.00 

21000 

0 

-.00 

1000 

2000 

103233.12 

3000 

2000 

-220.84 

5000 

2000 

-2687.20 

7000 

2000 

-1469.13 

9000 

2000 

-694.81 

11000 

2000 

-323.74 

13000 

2000 

-152.78 

15000 

2000 

-73.20 

17000 

2000 

-35.29 

19000 

2000 

-16.79 

21000 

2000 

-7.61 

1000 

4000 

-106052.11 

3000 

4000 

-22433 32 

5000 

4000 

-6881.11 

7000 

4000 

-2651.36 

9000 

4000 

-1133.96 

11000 

4000 

-514.44 

13000 

4000 

-242.18 

15000 

4000 

-116.43 

17000 

4000 

-56.17 

19000 

4000 

-26.49 

21000 

4000 

-11.62 

1000 

6000 

-112747.42 

3000 

6000 

-24307.53 

5000 

6000 

-6793.48 


CFNTER 2 

CHARGE 

CENTER 3 

COULOMBS 

HEIGHT 

COULOMBS 

-36.80 

2000.00 

12.70 


Y FIELD 

TOTAL 

THETA . 

-19578.95 

19578.95 

-1 . 5 7 08 

3637.27 

3637.27 

1.5708 

3716.05 

3716.05 

1.5708 

1 895. 86 

1895.86 

1.5708 

834.24 

834.24 

1.5708 

318.36 

318.36 

' 1 . 5708 

74.38 

' 74.38 

1.5708 

-41.76 

41.76 

-1.5708 

-97.60 

-97.60 

-1 5708 

-124.44 

-124.44 

-1.5708 

-137.05 

137.05 

-1.5708 

27386.64 

106804.06 

.2593 

11093.19 

11095.39 

1.5907 

4168.82 

4959.85 

2.1434 

1723.87 

2264.96 

2.2766 

750.14 

1022.49 

2.3179 

326.-33 

459.67 

2.3522 

131 .82 

201.79 

2.4297 

39.62 

83.24 

2.6456 

-4.78 

35.62 

-3.0069 

-26.07 

31.01 

-2.1430 

-35.92 

36.72 

-1.7795 

132575 97 

169774.69 

2.2455 

12080.09 

25479.06 

2.6476 

2186.94 

7220.28 

2.8339 

604.03 

2719.30 

2.9176 

219.37 

1154.98 

2.9505 

83.59 

521.19 

2.9805 

22.49 

243.22 

3.0490 

-7.75 

116.69 

-3.0752 

-22.75 

60.60 

-2.7568 

-29.72 

39.81 

-2.2988 

-32.40 

34.42 

-1.9151 

-120485.62 

1.65011 .41 

-2.3230 

-12275.94 

27231 .50 

-2.6739 

-3155.99 

7490.77 

-2.7067 
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• X 

Y 

X FIELD 

Y FIELD 

TOTAL 

THETA 

7000 

6000 

-2468.76 

-1113.56 

2708.28 

-2.7179 

9000 

6000 

-1045.96 

-455.83 

1140.97 

-2.7306 

n ooo 

6000 

-481.72 

-211.68 

526.18 

-2.7276 

13000 

6000 

-231.42 

-113.86 

257.91 

-2.6844 

15000 

6000 

-112.84 

-72.00 

133.85 

-2.5737 

17000 

6000 

-54.37 

-52.51 

75.59 

-2.3736 

19000 

6000 

-24.85 

-42.35 

49.10 

-2.1015 

21000 

6000 

-9.82 

-36.26 

37.57 

-1 .8354 _ 

1000 

8000 

6212.02 

-60396.07 

60714.70 

-1 .4683 

3000 

8000 

-268.39 

-18165.98 

18167.96 

-1.5856 

5000 

8000 

-1294.72 

-5611 .02 

5758.46 

-1.7976 

7000 

8000 

-840.03 

-2057.74 

2222.60 

-1.9584 

9000 

8000 

-468.22 

-867.77 

986.04 

-2.0656 

11000 

8000 

-251 .55 

-408.78 

479.98 

-2.1225 

13000 

8000 

-131.84 

-212.91 

250.42 

-2.1252 

15000 

8000 

-66.51 

-122.65 

319.53 

-2.0677 

17000 

8000 

-31 .09 

-78.18 

' 84.14 

-1.9493 

19000 

8000 

-12.11 

-54.75 

56.07 

-1.7885 

21000 

8000 

-2.16 

-41 .49 

41.54 

-1 .6227' 

1000 

10000 • 

178385.59 

-10509.62 

178694.91 

-.0588 

3000 

10000 

15647.43 

-6656.05 

17004.27 

-.4022 

5000 

10000 

3283.88 

-3374.60 

4708.69 

-.7990 

7000 

10000 

794.46 

-1638.40 

1802.85 

-1.1193 

9000 

10000 

188.63 

-817 80 

839.28 

-1 3441 

11000 

10000 

36.77 

-428.97 

430.54 

-1.4853 

13000 

10000 

3.73 

-238.40 

238.43 

-1 .5552 

15000 

10000 

1 .04 

-141 .21 

141 .22 

-1.5634 

17000 

10000 

4.47 

-89.55 

89.66 

-1.5209 

19000 

10000 

7.70 

-60.86 

61.35 

-1.4449 

21000 

10000 

9.56 

-44.18 

45.20 

-1.3577 

1000 

12000 

15343.80 

27077 87 

31123.04 

1 .0553 

3000 

12000 

9643.77 

3652.56 

10312.30 

.3621 

5000 

12000 

3639.16 

-331 55 

3654.24 

-.0909 

7000 

12000 

1391.71 

-660.15 

1540.34 

-.4429 

9000 

12000 

562.67 

-487.56 

744 52 

-.7140 

11000 

12000 

244.91 

-314.33 

398.48 

-.9089 

13000 

12000 

117.72 

-198.32 

230.62 

-1 .0351 

15000 

12000 

64.31 

-127.23 

142.56 

-1.1028 

17000 

12000 

40.48 

-84.50 

93.70 

-1.1241 

19000 

12000 

28.86 

-58.65 

65.36 

-1.1135 

21000 

12000 

22.48 

-42.67 

48.23 

-1.0860 

1000 

14000 

2204.57 

7239.22 

7567.46 

1.2752 

3000 

14000 

3376.22 

3169.39 

4630.75 

.7538 

5000 

14000 

2220.67 

792.93 

2357.99 

.3430 

7000 

14000 

1197.03 

16.01 

1197.14 

.0134 

9000 

14000 

620.33 

-159.92 

640.61 

-.2523 

11000 

14000 

327.34 

-162.03 

365.25 

-.4596 

13000 

14000 

181.07 

-127.25 

221.31 

-.6126 
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X 

Y 

X FIELD 

15000 

14000 

106.90 

17000 

14000 

67.96 

19000 

14000 

46.49 

21000 

14000 

33.92 

1000 

16000 

604.37 

3000 

16000' 

1259.94 

5000 

16000 

1160.41 

7000 

16000 

815.26 

9000 

16000 

516.19 

11000 

16000 

317.21 

13000 

16000 

196.48 

15000 

16000 

125.36 

17000 

16000 

83.34 

19000 

16000 

57.98 

21000 

16000 

42.17 


r FIELD 

TOTAL 

THETA 

-93.18 

141.81 

-.7170 

-67.41 

95.72 

-.7814 

-49.41 

67.85 

-.8158 

-37.15 

50.31 

-.8308 

2822.63 

2886.61 

1.3599 

1806.48 

2202.46 

.9618 

819.86 

1420.81 

.6151 

271.42 

859.25 

.3214 

39.06 

517.67 

.0755 

-40.16 

319.74 

-.1259 

-57.73 

204.79 

-.2858 

-54.10 

136.53 

-.4074 

-45.11 

94.76 

-.4961 

-36.24 

68.37 

-.5586 

-28.97 

51.16 

-.6010 




APPENDIX B 


The following program, written for the HP2100 computer, is used to 
minimize the percentage of .variation as shown in Equation 8 by 
varying the values of three charge centers A, B and C which are 
located at VA, YB and YC. One charge center is modified at a time 
beginning with A, then B, then C, back to A and so on until an 
oscillation is noticed in the percentage of variation. At that point, 
the value of INC can be decreased and the program restarted using the 
A, B and C values which were computed as the oscillation began. The 
output gives the values of A, B and C in coulombs and the corresponding 
sum of the squares of the percentage of variations. 

The initial value of INC was 0.1. Oscillation using this value of 
INC occurred at A = 6.8, B = 13.8 and C = 3.4 for the Simpson and 
Robinson model. ESQRT equalled .187 x 10 4 . Oscillation in the Gish 
and Wait model occurred at A = 2.7 and B = 5 with ESQRT equalling 

3 

.813 x 10 while in the Malan model it occurred at A = 22.9, 

B = -40, C = 10 and ESQRT = .269 x 10 4 . 

Changing INC to equal .01 and using the A, B and C values found in 
the previous step resulted in the following values for the modified 
models: 
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MODIFIED MODEL 

A • 

B 

C 

ESQRT 

Simpson & Robinson 

4.6 

-11.3 

5.97 

522 

Gish & Wait 

1.52 

- 2.64 

0.0 

624 

Mai an 

i 

20.1 

-36.8 

12.7 

799 
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PROGRAM FIELD 
REAL INC, MULT 

DIMENSION EFLD(9) ,X(9) ,EF(9) ,Y(9) 

11 READ (5,3) A,D,C,YA,YB,YC 

3 FORMAT (6F10.1 ) 

IF (A) 13,12,13 

13 READ (5,4) A1 ,A2,A3,A4,A5,A6,A7 

4 FORMAT (7A4) 

WRITE (6,5) A1 ,A2,A3,A4,A5,A6,A7 

5 FORMAT (1H1,45X,7A4,//',3X, "CHARGE CENTER A" ,5X, "CHARGE CENTER B", 5 
1 X , "CHARGE CENTER C " ,/ , 2X , "HEIGHT" , 3X , "COULOMBS " , 3X , "HEIGHT" , 3X , "CO 
2ULOMBS" , 3X , " HEIGHT" , 3X , "COULOMBS" ) 

WRITE (6,3) YA,A,YB,B,YC,C 

VALUE=A 

INC- J 

J=1 

EF ( 1 )=-1310. 

EF(2)=-1160. 

EF( 3 )=— 760 . 

EF(4)=-460. 

EF(5)=-260. 

EF(6)=— 1 35 . 

EF(7)=-60. 

EF(8)=1 . 

EF(9}=33000. 

X(1 )=2500 . 

X(2)=5000. 

X(3)=7500. 

X(4)=10000. 

X(5)=l 2500 . 

X(6)=15000. 

X( 7 ) =1 7500 . 

X (8) =20000. 

X( 9 ) = 1 31 00. 

MULT=1 ./(4.*3.141 5*8. 854* (10.** (-12))) 

131 11=1 

K=1 

130 EXQRT=0. 

DO 10 1=1 ,8,1 

EFLD(I) =MULT* ( ( 2 . *YA*A/ ( X ( I )**2 .+YA**2 . )**1 . 5 )+ ( 2. W YB*B/ ( X ( I )**2 .+ 
1YB**2. )**1 .5)+(2.~YC*C/(X(I)**2.+YC**2.)**l .5)) 

ESQRT=ESQRT+(EF(I )‘-EFLD(I ) )*(£F(I )-EFLD(I) )/(EF(I )*EF(I ) ) 

10 CONTINUE 

EFLD(9)=MULT*(A/ (Y(9)-YA)* W 2.-A/(Y(9)+YA)**2.+B/(Y(9)-YB)**2.-B/(Y 
1(9)+YB)**2.+C/(Y(9')-YC)**2.-C/(Y(9)+YC)**2.) 
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REPEODUOBILmOFfflB 


ESQRT=ESQRT+(EF(9)-EFLD(9) )*(EF(9)-EFL0(9) )/ (EF(9)*EF(9) ) 
ESQRT=ESQRT*100. 

IF (II-l) 20,20,30 
20 ESMIN=ESQRT 
GO TO 40 

30 I-F (ESQRT-ESMIN) 20,20,50 
50 IF ( I 1-2) 60,60,70 
60 K=-K 

VALUE=VALUE+K*INC 
40 11=11+1 

VALUE=VALUE+K*INC 
IF (ESMIN-100. ) 180,180,80 
80 IF (J-l) 90,90,100 
90 A=VALUE 
GO TO 130 

100 IF (0-2) 110,110,120 
110 B= VALUE 
GO TO 130 
120 C=VALUE 
GO TO 130 

70 IF (0-1) 140,140,150 
140 A=VALUE 
VALUE=B 
0=2 

GO TO 190 

150 IF (0-2) 160,160,170 
160 B=VALUE 
VALUE=C 
0=3 

GO TO 190 
170 C=VALUE 
VALUE=A 
0=1 

190 WRITE (6,2) A,B,C,ESMIN 
2 FORMAT (4E10.3) 

GO TO 131 

180 WRITE (6,6) , ^ 

6 F0RMAT(//, IX, "CHARGE A" ,2X, "CHARGE B" ,2X , "CHARGE C" ,3X ,"EFLD(1 ) " ,3 
1X,"EFLD(2)" ,3X,"EFLD(3)" ,3x,"EFLD( A )" ,3X ,"EFLD(5) " ,3X, U EFLD(6)" ,3X 
2, U EFLD(7)",3X,"EFL0(8)",3X,"EFL0(9)") 

WRITE (6,1 ) A,B,C,EFLD(1) ,EFLD(2) ,EFLD(3) ,EFL0(4) ,EFLD(5) ,£FLD(6) , 
1EFLD(7) ,EFLD(8) ,EFLD(9) 

1 FORMAT (12E10.3) 

GO TO 11 
12 END 
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INTRODUCTION 


OF 'itua 


SSSSEf® 1 ^ 


An alternative to passive sound insulation is active electronic control. 

In the approach, described earlier by Olson [1,2] and others [3,4], one 
connects a microphone, an inverting amplifier, and a loudspeaker in a 
negative feedback fashion to reduce noise. In this system (figure 1), 
the electric signal produced by the microphone from an incident sound 
pressure P , is amplified by an inverting amplifier with a gain of -K a - 
The amplifier drives the loudspeaker, producing an output pressure P . 

This pressure is summed acoustically with the incident sound pressure P. . 
The resulting pressure at the microphone is the sound pressure error 
signal, P £ ; 

= Pi + P„ 

£ 10 

The resulting sound pressure level P £ at the microphone has been shown to 
be less than the incident sound pressure level P. ] , [1,2,5]. If the 
microphone is placed near one's ear, a considerable reduction in the 
noise level is observed. However, airborne noise reduction through active 
electronic control is confined to a fairly narrow range of frequencies. 

This treatise will show that signal time delays from wave propagation 
through the acoustic medium and inherent transducer delay are responsible 
for this narrow bandwidth limitation. A classical method of lag- lead 
compensation will be used to improve the bandwidth of airborne-noise- 
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reduction systems developed in the references [1,2,5]. Finally, lag-lead 
compensation using maximally-flat (Butterworth) and equal-ripple (Chebyshev) 
polynomials will be developed to further extend the bandwidth. 

MATH MODEL 


The active compensation system of figure 1 is depicted schematically by 
the block diagram in figure 2. Here, -K represents the gain of an ideal 
amplifier. G m (s), G^(s), and G s (s) are the transfer functions of the 
microphone, filter, and loudspeaker respectively , in the s domain. 

(The transfer function G^(s) is the principle topic of this paper and 
is discussed below m detail.) System time delay x is represented by 
e ST . T' is transformation from the time domain to the frequency domain 
has been accomplished according to the Laplace shift theorem: 

L(y(t-T) u (t)} = e~ ST Y(s), , m 

where U (t) = 10 t<x 

T ( 1 t >x 

T “ T n + T t s where x = wave propagation delay, 

P P (3) 

and x^ = total transducer delay. 

Time delay x is the sum of transducer delay and acoustic wave propagation 
delay. When time delay is present in a feedback control system, numerical 
analysis of that system is difficult because of the introduction of the 
transcendental function e ST , m the transfer function. However, a 
graphical approach is often adequate for stability analysis of feedback 
systems with time delay. This is the method used here. 




Active compensation scheme 
Figure 1 
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The open-loop response (where G^(s) =1) of an active electronic 
airborne-noise controller are shown in figures 3a and 3b. At the higher 
frequencies, the phase-frequency relationship becomes linear, indicating 
the presence of time delay. This is graphically illustrated in figure 4. 
Note that the system phase is the sum of the phase associated with two 
poles and the phase that represents pure time delay. At high frequencies, 
for v/hich all phase associated with poles and zeros can be considered a 
constant, the delay t of equation (4) is simply the slope of the system 
phase function. 

4e“ JWX = -cor (4) 

For the system represented in figures 3a, 3b, and 4 a time dela> of 
49 psecs is determined. Meeker [5] achieved this result with a 
microphone-speaker separation of 3.1 mm (.125 inch), which accounts 
only for 9 psecs of signal delay due to acoustic wave propagation (with 
an acoustic wave velocity of 3^7 meters/sec.) Approximately 40 psecs 
of delay must be attributed to the transducers. This conclusion is 
substantiated by the work of C.A. Ewaskio and O.K. Mawardi [6]. Their 
measurements of phase shift in loudspeakers indicated time delays of the 
magnitude found above. Transducer delay in equation (3) is equal to 
40 psecs for the system of figures 3a and 3b. 



magnitude 
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The magnitude function in figure 3a indicates a bandwidth of 4 KH Z , at 
which point it rolls-off at -12 dB/octave. For this sytem to achieve 
stability without compensation, the open-loop gain must be made less than 
unity. 

Lowering the gain of the system, however, decreases its effectiveness. 

The acoustic noise intensity I is calculated by comparing the sound 
pressure error signal P £ with the original incident sound pressure P 1 
I in decibels is: 

PM<a) 

M J “> - 20 log 10 (6) 

The sound pressure is measured at the microphone ‘s diaphragm, P £ is 
measured with the active electronic control functioning, P is measured 
with the system off. 


From the block diagram (figure 2), P £ for real frequencies is. 


P„(j«) = 


P^jw) 


1+ !< a G t (jw)e‘ 


where G t (j«) = G m (jo>)G f ( jw)G s ( ju) . 


Substituting into equation (5), I becomes 
I r (jco) = 20 log-jQ | l+K a G t (j W ) e‘ jwT |. 


( 7 ) 
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For those real frequencies w, where the open-loop gain is much greater 
than unity, j K G t (jw) e~ JWT j »1 , the reduction in the acoustic noise 
intensity at the microphone's diaphragm is approximately equal to the 
system's open-loop gain. 

I r (jw) = 20 log-jQ |K a G t {jw) | , since [ e" Jt0T i =1. (8) 

Thus, for maximum noise reduction, the open-loop gain should be kept as 
high as possible. 

SYSTEM DESIGN 


To obtain both stability and high gain, bandwidth of the system must be 
sacrificed. This conclusion will be nrc.de clear through graphical stability 
analysis in the form of a Nyquist diagram. A polor plot of the open-loop 
frequency response of K G*( jw)e~’- ,il)T for is constructed. This plot 

at ■*-“ — * 

is referred to as a Nyquist diagram of K,G + .(jw)e -,:,a)T > see figure 3b. A 

at 

minimal-phase (no right half-plane zeros) closed-loop system is stable if 
and only if the Nyquist diagram of K,G^.( jcole - ^ 01 does not encircle the 
(-1 ,j0) critical point on the complex plane. 

Olson [1,2] and Meeker [5] used lag compensation to obtain stability, 
simultaneously obtaining high open-loop gain over a narrow range of 
frequencies, see figures 5a and 5b. To obtain a gain of 20dB, the 
bandwidth had' to be lowered to 200Hz. The compensated noise-reducing 
system of figure 4a has a filter G f (jw), where 


" jw/1250 +1 


(9) 
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The filter G f (jw) has one pole at 200 Hz. Thus, the Nyquist diagram 
avoids encirclement of the (-1 ,jO)critical point, producing a stable system. 

From equation (8) and figure 5a, reduction in the noise intensity is found 
to be approximately 20dB up to a frequency of 200 Hz. To increase the 
bandwidth of the system in figures 4a and 4b one would like zo raise the 
crossover frequency of ( w _-jgg ^ 5 defined to' be the 

frequency where Ae~^ Wt i's equal to -180 degrees). One would like to move 
w -180 1:0 that 1s eliminate the time delay; however, this is not 

physically realizable. In the airborne-noise-reduction system investigated, 
a large percentage of the time delay is associated with transducers. With 
present state-of-the-art transducers, little can be done to reduce the 
inherent time delay. Therefore, to avoid encirclement of the (-l,j0 ) 
point on the Nyquist diagram, the magnitude function must be less than 
unity before the frequency i qq * Consequently, this frequency, w_-jgg 
must be considered the upper limit on an airborne-noise-reduction system's 
bandwidth. 

If the bandwidth is to be increased, it must be done through manipulation 

of the magnitude function. An increase in the bandwidth may be realized 

if the magnitude function is made to roll-off at a greater rate than that 

of the one-pole filter in figure 5a. This allows the -3dB bandwidth to 

move closer to the upper limit, w_-jgQ. In a minimal-phase system 
increased roll-off is always accompanied by additional phase. Hence, 
trade-offs must be made. 
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Lag-lead compensation is commonly used to reduce gain while minimizing 
excess phase. The transfer function of a lag-lead compensator consists 
of a simple pole-zero pair: 

l+j“/co 7 

G f (jco) = 2- • (10) 

l+jto/Wp 

The pole is always nearer the origin of the s plane than the zero. Note 
that the magnitude function, G^(jw), has an initial gain of unity. As 
w approaches infinity, the gain becomes . Moreover, the phase 

function returns to zero as w goes to infinity. Lag- lead compensation 
allows the magnitude function to decrease while introducing smaller amounts 
of excess phase than lag compensation. 

The use of a first-order (both denominator and numerator are first order 
polynomials) lag-lead compensator in the airborne-noise-reduction systems 
of the references [1,2,5] would extend the bandwidth to 300Hz, see 
figures 6a and 6b. The pole-zero location of the compensator have been 
determined graphically and found to be 

l+jw/2.51x!0 4 

Gi(jco) = (11) 

I+ju/1880 

A graphical approach for finding the pole and zero locations of a lag-lead 
compensator is quite adequate for a single-input, single-output closed-loop 
system. The location of the pole and zero is left to the designer. The 




Nyquist diagram of Figure 6a 


Figure 6b 
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p-oper phase and gain margins may then be selected for desired overshoot 
and damping. Insofar as the noise-reduction system can be approximated 
by a second-order system, the phase-gain margin concept will be a 
grapm'cal convenience which can provide approximate closed-loop analysis 
by inspection. Criteria for phase and gain margins will not be 
discussed, except for the statement that phase margin greater than 
55 degrees and a gam margin greater than 6dB are desirable. This gives 
an approximate closed-loop damping >.55. 

The ideal normalized (the -3dB bandwidth frequency, defined as 
is equal to one) lag-lead compensator is characterized by a magnitude 
that is unity for frequencies, w<l and can be represented by 1/K^, for 
frequencies a.»1 . These ideal properties are not attainable with 
constant, lumped, and linear networks, therefore the requirements 
imposed by the ideal characteristics can be only approximated. This is 
done by allowing the magnitude and phase to stay within prescribed limits 
of the ideal, see figures 7a and 7b. 

It is evident from the magnitude function graphed in figures 8a and 8b 
that cascading more than one first-order lag-lead compensator results 
in a magnitude function [ G^(joj) | that tends away from the ideal. Also 
from figure 8b, it is evident that excess phase is increased with 
increased number of cascaded filters. 








phase magnitude 


OdB 



Figure 8a 



frequency 
Figure 8b 
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BUTTERWORTH-CHEBYSHEV LAG-LEAD COMPENSATION 


Functions of the form 


G f (s) 


N(s) 

W » 


02 ) 


which obey the specified constraints of figures 7a and 7b, represent 
approximations to the ideal lag- lead filter. Combinations of the 
maximal 1 y— f 1 at (Butterworth) and the equal-ripple (Chebyshev) polynomials 
provide candidate approximations to the ideal. 


The excess phase associated with the maximally- flat function increases 
with frequency at a much smaller rate than that of the equal -ripple 
polynomial, see figure 9aand9b. Of the four ratio combinations, 
representing the equal-ripple polynomial by N(s) (which contributes 
positive phase) and the maximal ly-flat polynomial by D(s) (which 
contributes relative little negative phase) yields the most desirable 
phase function. 


A suitable representation of this lag-lead compensator can be obtained 
by starting withthe general expression of the magnitude-squared function. 
Let 


|6f(j«) ! 2 = 


N(cu 2 ) 

D(</) 


(13) 


for real frequencies w. 



magnitude 


if 


jnaxjimal ly-flat 


ecjual -ripple 


frequency (radians) 

Fourth-order Butterworth and Chebyshev magnitude function 

Figure 9a 


equal -rij 
P 


maxii 

nally-flai 

phase 



frequency (radians) 

Fourch-order Betterworth and Chebyshev phase function 


Figure 9b 
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For frequencies uk! the magnitude-squared function of the normalized 

filter G^(jw) should be approximately unity. For frequencies w>l the 

2 

function should not depart apprecia^y from 1/K^. The filter design to 
be discussed in detail is 


ii 2 

g -f(Tw) I 


l+ e 2 C 2 ( J_ } 
n ' aw ' 


K 2 + e 2 2 


2IHTT 


(aw) 


2n 


where the C n (x) are called the Chebyshev polynomials. 


(14) 


(— L) = cos 2 (n cos"^ -J— ) 
aw aw 


2 -11 

cosh (n cosh ) 

aw 


(15a) 

(15b) 


o 

Although either form of C n (l/ow) can be employed, it is convenient to 
use equation (15a) for frequencies l/aw<1 and equation (15b) for 1/ciw>l. 

, I* 2 

It will presently be shown that jG^^w) | is a ratio of an equal -ripple 

2 2 
polynomial (N(w )) to a maximal ly-flat polynomial (D(w )) 

For real frequencies w, when w-*0, 


.. 2 
iG f (jw) | 

w+o 


2 r 2 
e C n 


2(n-l ) 


2 2 


(aw) 


2n 


(16a) 


because 


2 

cosh (n 


cosh 


•1 1)5 2 2(n-l)/J_j 


2n 


aw 


aw 


for 


1 

aw 


> 1 


(16b) 
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For real frequencies w, when 

„ 2 / 9 

| G f ( jw) i = j 1 /n f^r n odd 

f (l+£ 2 )/£ for n even, 
when n is the order of' the polynomial. 

This is true because 

cos 2 (ncos" ] = cos 2 (n |) . for «1 . 

For frequencies w>l/a, C n (l/aw) is a cosine-squared function. 
0< |C 2 (l/oo))| <1. 

0 

The magnitude characteristic of N(to ) is equal-ripple m nature 

■I 2 2 

numerator of iG^jw) j will swing between 1 and 1+s exactly 
for w>l/a. 

Note that 

| N(w 2 ) j = 1 when C 2 ( l/oo)) = 0. 

"l tt[c ~ 

This occurs when cc = i. sec (where k = 1,3,5, ). 

a 

Also , 

|N(o) 2 ) | = 1+c 2 when C 2 (-~) = l._ 

"j #jy 

This occurs when (o = — sec — (where k=0,l ,2, • • • • ) • 

a n 


(17a) 

07b) 

Therefore 

(18) • 

, the 
times , 


Because the parameter £ controls the amount of ripple for w>l/a , it is 
called the ripple factor. 
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All filter designs discussed will be normalized about the filter's 
-3dB frequency, w_ 3d g; thus = 1. The normalizing parameter a is 

the value for which 

2 

! G f( J “-3dB) 1 = 1 09) 

The value of a must be solved by a successive-approximation. The value 
of a as a function of K^/e for different n is shown in figure 10. If the 
approximation (16b) is valid for w=l , then 


C^(i) * 2 2 < n -’> — .for i> 1. 


(a) 


2n 


1 

a 


Substituting for C^(l/a) in equation (14), 


|G f (j“) 


.2 2 2(n-l) _J 


a 


2n 


w-1 


T2(n T) 
+ £ ^2F“ 


1 
2 ' 


Solving for the normalizing parameter a. 


a = ( y ) 1/n 2^ n " 1)/n < 1 

Y "f 


( 20 ) 


This result is displayed in figure 10, where it is compared with the 

calculated value of the normalizing parameter cc. Mote that equation (20) 

is valid for a less than unity, because cosh(x) in the Chebyshev 

polynomial may be approximated by e x /2. Therefore, equation (20) can be 

considered valid when 

Kf> ?n-l 
e 


(21) 
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Also of interest is the frequency w= 1/a . |G^(j“) j becomes 


2 

" ~ 1 4 £ ^ 




( 22 ) 


a 


However, if inequality (21) is valid, then from (22) 

2 

lG f (jw)l , 

/.i— L 


1+e 2 


(23) 




a 




2 

The denominator 0(s) becomes a better approximation to for “>!/<*. 

2 

Since the numerator N(s) for these frequencies swings between 1 and l+e 

ii 2 

(see equation (18)). the magnitude-squared function j G^(jw) | will swing 
1/K 2 and (l+s?)/K 2 , n times. This ripple effect begins at w=l/a which may 
be considered the corner frequency of the numerator. The above characteristics 
are diagrammatically illustrated in figure 11. 


Heretofore, just the magnitude function has been studied, because phase 
information cannot be obtained from magnitude-squared functions. To 
acquire the phase functions, the roots of the equal-ripple and maximal 1y- 
flat polynomials must be found. Expressing the magnitude-squared function 
as 

2 

| G f (j«) j = G f (j<o)G f (-jw) (24a) 


or 


2 

|G f (jo))| = G f (s)G f (-s) 


-s 


2 


to 




(24b) 



magnitude 
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G. (s) may be found. It should be remembered that although the notation 

is in terms of the imaginary variable jw, the variable which appears in 

o 

the function is the real variable w (for example see (14)). The desired 
generalization is made by substituting s/j foru. The remaining problem 
is to separate equation (2^-b) into its tv/o constituents, G^(s) and G^(-s). 
This is done by realizing that G^(s)G^(-s) contains the poles and zeros 
of G^(jco) and their mirror images with respect to the jw axis. Therefore, 
to assure a stable minimal -phase system one assigns all the left half-plane 
roots to Gf(ja). 

Solving the denominator of equation (14) by substituting s/j for w, 
D(s)D(-s) becomes 

p p p2( n-1 ) 

D(s)D(-,) = 4^—2^ ~ 0 . (25a) 

(-1 ) a S 


To simplify the computation, a new variable p is defined: p=l/s. 

Substituting into (25a), 


.2n 


- ("D n 


a 


2 \ 2 


2 , 

£ 2 


2Tn^TT 


(25b) 


Solving for the 2n roots, 



-ft/ i 2m-l x 

f e j 2 U+ n > 


, where m 1 ,2,3, • • • ,2n. 


(25c) 
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The n left half-plane roots, p. , may be found from (25c) to be 
n 1 

P k ~~x t-sirv 2 ^{2k-l) + j cos^(2k-l)] (25d) 

V e2 

, where k=l ,2, - • • ,n. 


Note that if the 


The roots p t , lie on a circle of radius —* n / a K f . 
normalizing parameter a is approximated by* 5 (20) the circle becomes.a unit 


circle. 


The function D(s) is found by replacing p with 1/s; 

D(s) = C 0 [D(p) | ]. (25e) 

S T 

Since D(s) is constructed from the roots , and the roots solve all polynomials 
that are multiples of D(p), the constant C Q is needed, The constant C Q is 
found by equating D(s) with the square root of equation (25a) for any to. 

Solving for the roots of the numerator of equation (14) is more difficult. 
Substituting s/j for , 

N(s)N(-s) = Ke 2 C 2 (j/ots) - 0. (26a) 

Again utilizing p=l/s for simplification, equation (26a) can be written 

1+A2 (f) = 0. 


(26b) 
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With equation (15a), 


C*; (^) = cos(n cos" 1 ^-) = + j . 

n ' a a ' — J e 


In solving for p, the complex variable w is introduced, 


w = u+jv = cos 


- 1 IP 


a 


Substituting and equating real and imaginary parts, 
cos nw = cos(nu+jnv) 

= cos(nu)cosh(nv) - jsin(nu)sinh(nv) 


+ jl 

— J e 


Hence, 


and 


cos(nu)cosh(nv) = 0 
sm(nu)s‘'nn(nv) = + l- 


Cosh(nv) cannot be equal to zero for real v; therefore, cos(nu) 
solutions for u are 

u m = (2m-l) , where m = 1,2, ,2n. 


For u ffl , sin(nu m ) is equal to one. Solving the remaining “equation 


m 


smh(nv) = + — , it follows that v = r-Uinh" 1 — . 


— e 


Rearranging equation (26d) 
JPn 


m 


a = cosw = cos(u m +jv) 


where m=l ,2, • • • ,2n. 


(26c) 

(26d) 

(26e) 

(26f) 

0. The 

(26g) 

of (26f) 
(26h) 

(261) 
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The n left half-plane roots are, 

P k = a(-sin(u k )sTnh(v)+jcos(u k )cosh(v)) , (26j) 

where k= 1 ,2, . , . ,n. 


The function N(s) is found by 

N(s) = C-jCN(p) | 1 ]. 

S= P 

The constant C-j is needed for the same reasons that C Q was used in 
equation (25e). 

To realize the benefit of the filter design, one must apply the above 
results to the noise-reduction systems in the references [1,2,5]. For 
comparison with [1,2,5] the amplifier gain in figure 2 will be 20dB 

I I l 

The 6^ (ju) used for this example will be a third-order filter, having 
a filter gain, K^,, of -17dB and a ripple of 1.5dB. To obtain the value 
of e, il+7 must be made 1.5dB greater than 1, that is 

. 20 log 1Q V Hs 2 = 1.5dB. (27) 

Here sis found to be .6423. 

For a filter gain of -17dB, from equation (17a) is 

-17d B = 20 log 10 jjl- (28) 

f 

where is found to be 7.079. 
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In this example, inequality (21) will be considered adequate; then from 
(20) a is found to be .7133. 

These four parameters (n,a,e,K^) provide a good idea of what the normalized 

1 I l 

magnitude function j (jw) | will look like; the -3dB bandwidth frequency 
is at a)=l; the corner frequency, where the ripple effect begins is at 
oj=l/a; magnitude function will swing between 1/1^= .1413 and V 1+c^/K^ = .1679 
exactly n=3 times. These characteristics are illustrated in figure 11. 


Starting with equation (25d) the roots of the D(s) are found to be 
p-j = -.500 + j .866 , ?2 = -1.00, - -.500 - j.866. 

Multipling the roots together to form D(p), 

D(p) = (p+1 .00)[(P+.500) 2 +.866 2 ] 

= p 3 H2.00p 2 +2.00p+l . 


Using (25e), 
D(p) 


s= 


1 , 2.00 2.00 x n 
+ + _+ 1 

s s 


(29) 


Before. finding D(s), C of equation (25e) must be determined. As or*=° 

for real frequencies D(p)|^ - l(see (29)) but, |D(s)|^-K f (from the square 

s=l/p r 

foot of (25a)). Therefore, C Q must be equal to (see (25e)). It 
follows that 

D(s) = 7 .079[— L + + ZjQQ- + 1]. 


(30) 
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Using equation (26j) the roots are found to be 

p 1 = -.1495+0.6696, P £ = -.2990, P 3 =-.1495-0 .6696. 

Multipling the roots together to form N(p), 

N(p) = (p+.2990)[(p+.1^95) 2 +.6696 2 ] 

= p 3 +.5980p 2 +. 5601 p+. 1407. 




(32) 


(33) 


in 

The magnitude and phase functions of (joo) for real frequencies w are 
shown -in figure 12. 


II! 

G (s) is now used in the noise-reducing system in figures 3a and 3b 
f 

(with amplifer gain K a =20dB) to reduce the open-loop gain from to 
K /K*. Because the phase of G- (s) (for large n) can shift beyond -180 

a T T 

degrees for some frequencies, the open-loop gain |K a G t | must be kept 
greater than unity until the open-loop phase returns to above -180 
degrees. This is necessary to avoid encirclment of the (-1,j0) critical 



frequency (radians) 

Third-order Butterworth-Chebyshev Filter 
Figure 12 
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point on the Nyqinst diagram. Subsequently, the magnitude must become 
less than unity before the open-loop phase function, which reflects the 
increasing negative phase of the time delay, returns below -180 degrees, 
see 13a and 13b. 

This magnitude characteristic is not alv/ays guaranteed by G m (s)G s (s). 

In this case, an additional lag pole will be necessary. 

Through graphical manipulation (using figures 3a, 3b, 12a, 12b) , the 
open-loop transfer function characteristics in figure 13a and 13b were 

III 

constructed. Here, the -3dB frequency aj _ 3 C jg G^(jw) is found to be 
l.OKHz. The normalized filter function (33) is translated to the 
actual frequencies by replacing s in (33) by s/2fr(1000) : 

g"(s)| = 4.030xl0~ 12 s 3 -H.01xl0~ 8 s-i-6 76x1 0~ 4 s+7. 079 (M) 

2.85xl0~ 12 s 3 -i-3.59xl0 _7 s 2 +2 25x1 0‘ 3 s+7 .079 
w -3dB = 2000ir 

Also, from figures 13a and 13b, it can be seen that a pole is needed near 
2.2KHz to bring the magnitude function below unity. Now the open-loop 
transfer function in s for figures 13a and 13b is 

] 1 1 (35) 

7.23xl0~ 5 s+l (3.98x10' 5 s+1 ) 2 

-3db=20007r 

In the above example the bandv/idth was extended to 1000Hz with a third- 
order, maximally- flat, equal-ripple, lag-lead filter. 
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CONCLUSION 


The work reported upon here showed that, by using additional electronic 
filtering, the bandwidth of an active electronic noise reduction system 
could be increased by two octaves. This improvement was achieved 
through the use of the algebraical ly tractable Butterworth and Chebyshev 
functions. Also, filter parameters (n,o.,e,K^) were chosen such that 
significant filter characteristics could be determined by inspection. 

Hence, by employing such an easily designed electronic filter, a significant 
improvement in active electronic noise reduction systems was realized. 

For further improvements in system bandwidth, higher-order filter functions 
could be tried, in addition to types other than Butterworth and Chebyshev. 
Moreover, improvements which would reduce time delay in other system 
components should be investigated. 
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